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

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

 48:   Logically Collective

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

 55:   Level: advanced

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

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

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

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

 82: /*@
 83:   VecStrideNorm - Computes the norm of subvector of a vector defined
 84:   by a starting point and a stride.

 86:   Collective

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

 93:   Output Parameter:
 94: . nrm - the norm

 96:   Level: advanced

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

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

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

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

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

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

148: /*@
149:   VecStrideMax - Computes the maximum of subvector of a vector defined
150:   by a starting point and a stride and optionally its location.

152:   Collective

154:   Input Parameters:
155: + v     - the vector
156: - start - starting point of the subvector (defined by a stride)

158:   Output Parameters:
159: + idex - the location where the maximum occurred  (pass `NULL` if not required)
160: - nrm  - the maximum value in the subvector

162:   Level: advanced

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

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

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

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

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

211:     PetscCall(VecGetOwnershipRange(v, &rstart, NULL));
212:     in.v = max;
213:     in.i = rstart + id;
214:     PetscCall(MPIU_Allreduce(&in, &out, 1, MPIU_REAL_INT, MPIU_MAXLOC, PetscObjectComm((PetscObject)v)));
215:     *nrm  = out.v;
216:     *idex = out.i;
217:   }
218: #endif
219:   PetscFunctionReturn(PETSC_SUCCESS);
220: }

222: /*@
223:   VecStrideMin - Computes the minimum of subvector of a vector defined
224:   by a starting point and a stride and optionally its location.

226:   Collective

228:   Input Parameters:
229: + v     - the vector
230: - start - starting point of the subvector (defined by a stride)

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

236:   Level: advanced

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

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

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

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

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

285:     PetscCall(VecGetOwnershipRange(v, &rstart, NULL));
286:     in.v = min;
287:     in.i = rstart + id;
288:     PetscCall(MPIU_Allreduce(&in, &out, 1, MPIU_REAL_INT, MPIU_MINLOC, PetscObjectComm((PetscObject)v)));
289:     *nrm  = out.v;
290:     *idex = out.i;
291:   }
292: #endif
293:   PetscFunctionReturn(PETSC_SUCCESS);
294: }

296: /*@
297:   VecStrideSum - Computes the sum of subvector of a vector defined
298:   by a starting point and a stride.

300:   Collective

302:   Input Parameters:
303: + v     - the vector
304: - start - starting point of the subvector (defined by a stride)

306:   Output Parameter:
307: . sum - the sum

309:   Level: advanced

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

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

318: .seealso: `Vec`, `VecSum()`, `VecStrideGather()`, `VecStrideScatter()`, `VecStrideMin()`, `VecStrideMax()`
319: @*/
320: PetscErrorCode VecStrideSum(Vec v, PetscInt start, PetscScalar *sum)
321: {
322:   PetscInt           i, n, bs;
323:   const PetscScalar *x;
324:   PetscScalar        local_sum = 0.0;

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

340: /*@
341:   VecStrideScaleAll - Scales the subvectors of a vector defined
342:   by a starting point and a stride.

344:   Logically Collective

346:   Input Parameters:
347: + v      - the vector
348: - scales - values to multiply each subvector entry by

350:   Level: advanced

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

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

358: .seealso: `Vec`, `VecNorm()`, `VecStrideScale()`, `VecScale()`, `VecStrideGather()`, `VecStrideScatter()`, `VecStrideMin()`, `VecStrideMax()`
359: @*/
360: PetscErrorCode VecStrideScaleAll(Vec v, const PetscScalar *scales)
361: {
362:   PetscInt     i, j, n, bs;
363:   PetscScalar *x;

365:   PetscFunctionBegin;
367:   PetscAssertPointer(scales, 2);
368:   PetscCall(VecGetLocalSize(v, &n));
369:   PetscCall(VecGetBlockSize(v, &bs));
370:   PetscCall(VecGetArray(v, &x));
371:   /* need to provide optimized code for each bs */
372:   for (i = 0; i < n; i += bs) {
373:     for (j = 0; j < bs; j++) x[i + j] *= scales[j];
374:   }
375:   PetscCall(VecRestoreArray(v, &x));
376:   PetscFunctionReturn(PETSC_SUCCESS);
377: }

379: /*@
380:   VecStrideNormAll - Computes the norms of subvectors of a vector defined
381:   by a starting point and a stride.

383:   Collective

385:   Input Parameters:
386: + v     - the vector
387: - ntype - type of norm, one of `NORM_1`, `NORM_2`, `NORM_INFINITY`

389:   Output Parameter:
390: . nrm - the norms

392:   Level: advanced

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

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

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

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

405: .seealso: `Vec`, `VecNorm()`, `VecStrideGather()`, `VecStrideScatter()`, `VecStrideMin()`, `VecStrideMax()`
406: @*/
407: PetscErrorCode VecStrideNormAll(Vec v, NormType ntype, PetscReal nrm[])
408: {
409:   PetscInt           i, j, n, bs;
410:   const PetscScalar *x;
411:   PetscReal          tnorm[128];
412:   MPI_Comm           comm;

414:   PetscFunctionBegin;
417:   PetscAssertPointer(nrm, 3);
418:   PetscCall(VecGetLocalSize(v, &n));
419:   PetscCall(VecGetArrayRead(v, &x));
420:   PetscCall(PetscObjectGetComm((PetscObject)v, &comm));

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

425:   if (ntype == NORM_2) {
426:     PetscScalar sum[128];
427:     for (j = 0; j < bs; j++) sum[j] = 0.0;
428:     for (i = 0; i < n; i += bs) {
429:       for (j = 0; j < bs; j++) sum[j] += x[i + j] * (PetscConj(x[i + j]));
430:     }
431:     for (j = 0; j < bs; j++) tnorm[j] = PetscRealPart(sum[j]);

433:     PetscCall(MPIU_Allreduce(tnorm, nrm, bs, MPIU_REAL, MPIU_SUM, comm));
434:     for (j = 0; j < bs; j++) nrm[j] = PetscSqrtReal(nrm[j]);
435:   } else if (ntype == NORM_1) {
436:     for (j = 0; j < bs; j++) tnorm[j] = 0.0;

438:     for (i = 0; i < n; i += bs) {
439:       for (j = 0; j < bs; j++) tnorm[j] += PetscAbsScalar(x[i + j]);
440:     }

442:     PetscCall(MPIU_Allreduce(tnorm, nrm, bs, MPIU_REAL, MPIU_SUM, comm));
443:   } else if (ntype == NORM_INFINITY) {
444:     PetscReal tmp;
445:     for (j = 0; j < bs; j++) tnorm[j] = 0.0;

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

463: /*@
464:   VecStrideMaxAll - Computes the maximums of subvectors of a vector defined
465:   by a starting point and a stride and optionally its location.

467:   Collective

469:   Input Parameter:
470: . v - the vector

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

477:   Level: advanced

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

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

485: .seealso: `Vec`, `VecMax()`, `VecStrideNorm()`, `VecStrideGather()`, `VecStrideScatter()`, `VecStrideMin()`
486: @*/
487: PetscErrorCode VecStrideMaxAll(Vec v, PetscInt idex[], PetscReal nrm[])
488: {
489:   PetscInt           i, j, n, bs;
490:   const PetscScalar *x;
491:   PetscReal          max[128], tmp;
492:   MPI_Comm           comm;

494:   PetscFunctionBegin;
496:   PetscAssertPointer(nrm, 3);
497:   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");
498:   PetscCall(VecGetLocalSize(v, &n));
499:   PetscCall(VecGetArrayRead(v, &x));
500:   PetscCall(PetscObjectGetComm((PetscObject)v, &comm));

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

505:   if (!n) {
506:     for (j = 0; j < bs; j++) max[j] = PETSC_MIN_REAL;
507:   } else {
508:     for (j = 0; j < bs; j++) max[j] = PetscRealPart(x[j]);

510:     for (i = bs; i < n; i += bs) {
511:       for (j = 0; j < bs; j++) {
512:         if ((tmp = PetscRealPart(x[i + j])) > max[j]) max[j] = tmp;
513:       }
514:     }
515:   }
516:   PetscCall(MPIU_Allreduce(max, nrm, bs, MPIU_REAL, MPIU_MAX, comm));

518:   PetscCall(VecRestoreArrayRead(v, &x));
519:   PetscFunctionReturn(PETSC_SUCCESS);
520: }

522: /*@
523:   VecStrideMinAll - Computes the minimum of subvector of a vector defined
524:   by a starting point and a stride and optionally its location.

526:   Collective

528:   Input Parameter:
529: . v - the vector

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

536:   Level: advanced

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

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

544: .seealso: `Vec`, `VecMin()`, `VecStrideNorm()`, `VecStrideGather()`, `VecStrideScatter()`, `VecStrideMax()`
545: @*/
546: PetscErrorCode VecStrideMinAll(Vec v, PetscInt idex[], PetscReal nrm[])
547: {
548:   PetscInt           i, n, bs, j;
549:   const PetscScalar *x;
550:   PetscReal          min[128], tmp;
551:   MPI_Comm           comm;

553:   PetscFunctionBegin;
555:   PetscAssertPointer(nrm, 3);
556:   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");
557:   PetscCall(VecGetLocalSize(v, &n));
558:   PetscCall(VecGetArrayRead(v, &x));
559:   PetscCall(PetscObjectGetComm((PetscObject)v, &comm));

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

564:   if (!n) {
565:     for (j = 0; j < bs; j++) min[j] = PETSC_MAX_REAL;
566:   } else {
567:     for (j = 0; j < bs; j++) min[j] = PetscRealPart(x[j]);

569:     for (i = bs; i < n; i += bs) {
570:       for (j = 0; j < bs; j++) {
571:         if ((tmp = PetscRealPart(x[i + j])) < min[j]) min[j] = tmp;
572:       }
573:     }
574:   }
575:   PetscCall(MPIU_Allreduce(min, nrm, bs, MPIU_REAL, MPIU_MIN, comm));

577:   PetscCall(VecRestoreArrayRead(v, &x));
578:   PetscFunctionReturn(PETSC_SUCCESS);
579: }

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

584:   Collective

586:   Input Parameter:
587: . v - the vector

589:   Output Parameter:
590: . sums - the sums

592:   Level: advanced

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

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

601: .seealso: `Vec`, `VecSum()`, `VecStrideGather()`, `VecStrideScatter()`, `VecStrideMin()`, `VecStrideMax()`
602: @*/
603: PetscErrorCode VecStrideSumAll(Vec v, PetscScalar sums[])
604: {
605:   PetscInt           i, j, n, bs;
606:   const PetscScalar *x;
607:   PetscScalar        local_sums[128];
608:   MPI_Comm           comm;

610:   PetscFunctionBegin;
612:   PetscAssertPointer(sums, 2);
613:   PetscCall(VecGetLocalSize(v, &n));
614:   PetscCall(VecGetArrayRead(v, &x));
615:   PetscCall(PetscObjectGetComm((PetscObject)v, &comm));

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

620:   for (j = 0; j < bs; j++) local_sums[j] = 0.0;
621:   for (i = 0; i < n; i += bs) {
622:     for (j = 0; j < bs; j++) local_sums[j] += x[i + j];
623:   }
624:   PetscCall(MPIU_Allreduce(local_sums, sums, bs, MPIU_SCALAR, MPIU_SUM, comm));

626:   PetscCall(VecRestoreArrayRead(v, &x));
627:   PetscFunctionReturn(PETSC_SUCCESS);
628: }

630: /*----------------------------------------------------------------------------------------------*/
631: /*@
632:   VecStrideGatherAll - Gathers all the single components from a multi-component vector into
633:   separate vectors.

635:   Collective

637:   Input Parameters:
638: + v    - the vector
639: - addv - one of `ADD_VALUES`, `INSERT_VALUES`, `MAX_VALUES`

641:   Output Parameter:
642: . s - the location where the subvectors are stored

644:   Level: advanced

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

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

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

657:   Not optimized; could be easily

659: .seealso: `Vec`, `VecStrideNorm()`, `VecStrideScatter()`, `VecStrideMin()`, `VecStrideMax()`, `VecStrideGather()`,
660:           `VecStrideScatterAll()`
661: @*/
662: PetscErrorCode VecStrideGatherAll(Vec v, Vec s[], InsertMode addv)
663: {
664:   PetscInt           i, n, n2, bs, j, k, *bss = NULL, nv, jj, nvc;
665:   PetscScalar      **y;
666:   const PetscScalar *x;

668:   PetscFunctionBegin;
670:   PetscAssertPointer(s, 2);
672:   PetscCall(VecGetLocalSize(v, &n));
673:   PetscCall(VecGetLocalSize(s[0], &n2));
674:   PetscCall(VecGetArrayRead(v, &x));
675:   PetscCall(VecGetBlockSize(v, &bs));
676:   PetscCheck(bs > 0, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Input vector does not have a valid blocksize set");

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

691:   n = n / bs;

693:   jj = 0;
694:   if (addv == INSERT_VALUES) {
695:     for (j = 0; j < nv; j++) {
696:       for (k = 0; k < bss[j]; k++) {
697:         for (i = 0; i < n; i++) y[j][i * bss[j] + k] = x[bs * i + jj + k];
698:       }
699:       jj += bss[j];
700:     }
701:   } else if (addv == ADD_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: #if !defined(PETSC_USE_COMPLEX)
709:   } else if (addv == MAX_VALUES) {
710:     for (j = 0; j < nv; j++) {
711:       for (k = 0; k < bss[j]; k++) {
712:         for (i = 0; i < n; i++) y[j][i * bss[j] + k] = PetscMax(y[j][i * bss[j] + k], x[bs * i + jj + k]);
713:       }
714:       jj += bss[j];
715:     }
716: #endif
717:   } else SETERRQ(PetscObjectComm((PetscObject)v), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown norm type");

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

722:   PetscCall(PetscFree2(y, bss));
723:   PetscFunctionReturn(PETSC_SUCCESS);
724: }

726: /*@
727:   VecStrideScatterAll - Scatters all the single components from separate vectors into
728:   a multi-component vector.

730:   Collective

732:   Input Parameters:
733: + s    - the location where the subvectors are stored
734: - addv - one of `ADD_VALUES`, `INSERT_VALUES`, `MAX_VALUES`

736:   Output Parameter:
737: . v - the multicomponent vector

739:   Level: advanced

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

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

748:   Not optimized; could be easily

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

752: @*/
753: PetscErrorCode VecStrideScatterAll(Vec s[], Vec v, InsertMode addv)
754: {
755:   PetscInt            i, n, n2, bs, j, jj, k, *bss = NULL, nv, nvc;
756:   PetscScalar        *x;
757:   PetscScalar const **y;

759:   PetscFunctionBegin;
761:   PetscAssertPointer(s, 1);
763:   PetscCall(VecGetLocalSize(v, &n));
764:   PetscCall(VecGetLocalSize(s[0], &n2));
765:   PetscCall(VecGetArray(v, &x));
766:   PetscCall(VecGetBlockSize(v, &bs));
767:   PetscCheck(bs > 0, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Input vector does not have a valid blocksize set");

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

782:   n = n / bs;

784:   jj = 0;
785:   if (addv == INSERT_VALUES) {
786:     for (j = 0; j < nv; j++) {
787:       for (k = 0; k < bss[j]; k++) {
788:         for (i = 0; i < n; i++) x[bs * i + jj + k] = y[j][i * bss[j] + k];
789:       }
790:       jj += bss[j];
791:     }
792:   } else if (addv == ADD_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: #if !defined(PETSC_USE_COMPLEX)
800:   } else if (addv == MAX_VALUES) {
801:     for (j = 0; j < nv; j++) {
802:       for (k = 0; k < bss[j]; k++) {
803:         for (i = 0; i < n; i++) x[bs * i + jj + k] = PetscMax(x[bs * i + jj + k], y[j][i * bss[j] + k]);
804:       }
805:       jj += bss[j];
806:     }
807: #endif
808:   } else SETERRQ(PetscObjectComm((PetscObject)v), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown norm type");

810:   PetscCall(VecRestoreArray(v, &x));
811:   for (i = 0; i < nv; i++) PetscCall(VecRestoreArrayRead(s[i], &y[i]));
812:   PetscCall(PetscFree2(*(PetscScalar ***)&y, bss));
813:   PetscFunctionReturn(PETSC_SUCCESS);
814: }

816: /*@
817:   VecStrideGather - Gathers a single component from a multi-component vector into
818:   another vector.

820:   Collective

822:   Input Parameters:
823: + v     - the vector
824: . start - starting point of the subvector (defined by a stride)
825: - addv  - one of `ADD_VALUES`, `INSERT_VALUES`, `MAX_VALUES`

827:   Output Parameter:
828: . s - the location where the subvector is stored

830:   Level: advanced

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

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

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

842:   Not optimized; could be easily

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

859: /*@
860:   VecStrideScatter - Scatters a single component from a vector into a multi-component vector.

862:   Collective

864:   Input Parameters:
865: + s     - the single-component vector
866: . start - starting point of the subvector (defined by a stride)
867: - addv  - one of `ADD_VALUES`, `INSERT_VALUES`, `MAX_VALUES`

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

872:   Level: advanced

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

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

881:   Not optimized; could be easily

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

898: /*@
899:   VecStrideSubSetGather - Gathers a subset of components from a multi-component vector into
900:   another vector.

902:   Collective

904:   Input Parameters:
905: + v    - the vector
906: . nidx - the number of indices
907: . idxv - the indices of the components 0 <= idxv[0] ...idxv[nidx-1] < bs(v), they need not be sorted
908: . 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`
909: - addv - one of `ADD_VALUES`, `INSERT_VALUES`, `MAX_VALUES`

911:   Output Parameter:
912: . s - the location where the subvector is stored

914:   Level: advanced

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

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

922:   Not optimized; could be easily

924: .seealso: `Vec`, `VecStrideNorm()`, `VecStrideScatter()`, `VecStrideGather()`, `VecStrideSubSetScatter()`, `VecStrideMin()`, `VecStrideMax()`, `VecStrideGatherAll()`,
925:           `VecStrideScatterAll()`
926: @*/
927: PetscErrorCode VecStrideSubSetGather(Vec v, PetscInt nidx, const PetscInt idxv[], const PetscInt idxs[], Vec s, InsertMode addv)
928: {
929:   PetscFunctionBegin;
932:   if (nidx == PETSC_DETERMINE) nidx = s->map->bs;
933:   PetscUseTypeMethod(v, stridesubsetgather, nidx, idxv, idxs, s, addv);
934:   PetscFunctionReturn(PETSC_SUCCESS);
935: }

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

940:   Collective

942:   Input Parameters:
943: + s    - the smaller-component vector
944: . nidx - the number of indices in idx
945: . 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`
946: . idxv - the indices of the components in the larger-component vector, 0 <= idx[0] ...idx[nidx-1] < bs(v) they need not be sorted
947: - addv - one of `ADD_VALUES`, `INSERT_VALUES`, `MAX_VALUES`

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

952:   Level: advanced

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

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

960:   Not optimized; could be easily

962: .seealso: `Vec`, `VecStrideNorm()`, `VecStrideGather()`, `VecStrideSubSetGather()`, `VecStrideMin()`, `VecStrideMax()`, `VecStrideGatherAll()`,
963:           `VecStrideScatterAll()`
964: @*/
965: PetscErrorCode VecStrideSubSetScatter(Vec s, PetscInt nidx, const PetscInt idxs[], const PetscInt idxv[], Vec v, InsertMode addv)
966: {
967:   PetscFunctionBegin;
970:   if (nidx == PETSC_DETERMINE) nidx = s->map->bs;
971:   PetscCall((*v->ops->stridesubsetscatter)(s, nidx, idxs, idxv, v, addv));
972:   PetscFunctionReturn(PETSC_SUCCESS);
973: }

975: PetscErrorCode VecStrideGather_Default(Vec v, PetscInt start, Vec s, InsertMode addv)
976: {
977:   PetscInt           i, n, bs, ns;
978:   const PetscScalar *x;
979:   PetscScalar       *y;

981:   PetscFunctionBegin;
982:   PetscCall(VecGetLocalSize(v, &n));
983:   PetscCall(VecGetLocalSize(s, &ns));
984:   PetscCall(VecGetArrayRead(v, &x));
985:   PetscCall(VecGetArray(s, &y));

987:   bs = v->map->bs;
988:   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);
989:   x += start;
990:   n = n / bs;

992:   if (addv == INSERT_VALUES) {
993:     for (i = 0; i < n; i++) y[i] = x[bs * i];
994:   } else if (addv == ADD_VALUES) {
995:     for (i = 0; i < n; i++) y[i] += x[bs * i];
996: #if !defined(PETSC_USE_COMPLEX)
997:   } else if (addv == MAX_VALUES) {
998:     for (i = 0; i < n; i++) y[i] = PetscMax(y[i], x[bs * i]);
999: #endif
1000:   } else SETERRQ(PetscObjectComm((PetscObject)v), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown norm type");

1002:   PetscCall(VecRestoreArrayRead(v, &x));
1003:   PetscCall(VecRestoreArray(s, &y));
1004:   PetscFunctionReturn(PETSC_SUCCESS);
1005: }

1007: PetscErrorCode VecStrideScatter_Default(Vec s, PetscInt start, Vec v, InsertMode addv)
1008: {
1009:   PetscInt           i, n, bs, ns;
1010:   PetscScalar       *x;
1011:   const PetscScalar *y;

1013:   PetscFunctionBegin;
1014:   PetscCall(VecGetLocalSize(v, &n));
1015:   PetscCall(VecGetLocalSize(s, &ns));
1016:   PetscCall(VecGetArray(v, &x));
1017:   PetscCall(VecGetArrayRead(s, &y));

1019:   bs = v->map->bs;
1020:   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);
1021:   x += start;
1022:   n = n / bs;

1024:   if (addv == INSERT_VALUES) {
1025:     for (i = 0; i < n; i++) x[bs * i] = y[i];
1026:   } else if (addv == ADD_VALUES) {
1027:     for (i = 0; i < n; i++) x[bs * i] += y[i];
1028: #if !defined(PETSC_USE_COMPLEX)
1029:   } else if (addv == MAX_VALUES) {
1030:     for (i = 0; i < n; i++) x[bs * i] = PetscMax(y[i], x[bs * i]);
1031: #endif
1032:   } else SETERRQ(PetscObjectComm((PetscObject)s), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown norm type");

1034:   PetscCall(VecRestoreArray(v, &x));
1035:   PetscCall(VecRestoreArrayRead(s, &y));
1036:   PetscFunctionReturn(PETSC_SUCCESS);
1037: }

1039: PetscErrorCode VecStrideSubSetGather_Default(Vec v, PetscInt nidx, const PetscInt idxv[], const PetscInt idxs[], Vec s, InsertMode addv)
1040: {
1041:   PetscInt           i, j, n, bs, bss, ns;
1042:   const PetscScalar *x;
1043:   PetscScalar       *y;

1045:   PetscFunctionBegin;
1046:   PetscCall(VecGetLocalSize(v, &n));
1047:   PetscCall(VecGetLocalSize(s, &ns));
1048:   PetscCall(VecGetArrayRead(v, &x));
1049:   PetscCall(VecGetArray(s, &y));

1051:   bs  = v->map->bs;
1052:   bss = s->map->bs;
1053:   n   = n / bs;

1055:   if (PetscDefined(USE_DEBUG)) {
1056:     PetscCheck(n == ns / bss, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Incompatible layout of vectors");
1057:     for (j = 0; j < nidx; j++) {
1058:       PetscCheck(idxv[j] >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "idx[%" PetscInt_FMT "] %" PetscInt_FMT " is negative", j, idxv[j]);
1059:       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);
1060:     }
1061:     PetscCheck(idxs || bss == nidx, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Must provide idxs when not gathering into all locations");
1062:   }

1064:   if (addv == INSERT_VALUES) {
1065:     if (!idxs) {
1066:       for (i = 0; i < n; i++) {
1067:         for (j = 0; j < bss; j++) y[bss * i + j] = x[bs * i + idxv[j]];
1068:       }
1069:     } else {
1070:       for (i = 0; i < n; i++) {
1071:         for (j = 0; j < bss; j++) y[bss * i + idxs[j]] = x[bs * i + idxv[j]];
1072:       }
1073:     }
1074:   } else if (addv == ADD_VALUES) {
1075:     if (!idxs) {
1076:       for (i = 0; i < n; i++) {
1077:         for (j = 0; j < bss; j++) y[bss * i + j] += x[bs * i + idxv[j]];
1078:       }
1079:     } else {
1080:       for (i = 0; i < n; i++) {
1081:         for (j = 0; j < bss; j++) y[bss * i + idxs[j]] += x[bs * i + idxv[j]];
1082:       }
1083:     }
1084: #if !defined(PETSC_USE_COMPLEX)
1085:   } else if (addv == MAX_VALUES) {
1086:     if (!idxs) {
1087:       for (i = 0; i < n; i++) {
1088:         for (j = 0; j < bss; j++) y[bss * i + j] = PetscMax(y[bss * i + j], x[bs * i + idxv[j]]);
1089:       }
1090:     } else {
1091:       for (i = 0; i < n; i++) {
1092:         for (j = 0; j < bss; j++) y[bss * i + idxs[j]] = PetscMax(y[bss * i + idxs[j]], x[bs * i + idxv[j]]);
1093:       }
1094:     }
1095: #endif
1096:   } else SETERRQ(PetscObjectComm((PetscObject)v), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown norm type");

1098:   PetscCall(VecRestoreArrayRead(v, &x));
1099:   PetscCall(VecRestoreArray(s, &y));
1100:   PetscFunctionReturn(PETSC_SUCCESS);
1101: }

1103: PetscErrorCode VecStrideSubSetScatter_Default(Vec s, PetscInt nidx, const PetscInt idxs[], const PetscInt idxv[], Vec v, InsertMode addv)
1104: {
1105:   PetscInt           j, i, n, bs, ns, bss;
1106:   PetscScalar       *x;
1107:   const PetscScalar *y;

1109:   PetscFunctionBegin;
1110:   PetscCall(VecGetLocalSize(v, &n));
1111:   PetscCall(VecGetLocalSize(s, &ns));
1112:   PetscCall(VecGetArray(v, &x));
1113:   PetscCall(VecGetArrayRead(s, &y));

1115:   bs  = v->map->bs;
1116:   bss = s->map->bs;
1117:   n   = n / bs;

1119:   if (PetscDefined(USE_DEBUG)) {
1120:     PetscCheck(n == ns / bss, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Incompatible layout of vectors");
1121:     for (j = 0; j < bss; j++) {
1122:       if (idxs) {
1123:         PetscCheck(idxs[j] >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "idx[%" PetscInt_FMT "] %" PetscInt_FMT " is negative", j, idxs[j]);
1124:         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);
1125:       }
1126:     }
1127:     PetscCheck(idxs || bss == nidx, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Must provide idxs when not scattering from all locations");
1128:   }

1130:   if (addv == INSERT_VALUES) {
1131:     if (!idxs) {
1132:       for (i = 0; i < n; i++) {
1133:         for (j = 0; j < bss; j++) x[bs * i + idxv[j]] = y[bss * i + j];
1134:       }
1135:     } else {
1136:       for (i = 0; i < n; i++) {
1137:         for (j = 0; j < bss; j++) x[bs * i + idxv[j]] = y[bss * i + idxs[j]];
1138:       }
1139:     }
1140:   } else if (addv == ADD_VALUES) {
1141:     if (!idxs) {
1142:       for (i = 0; i < n; i++) {
1143:         for (j = 0; j < bss; j++) x[bs * i + idxv[j]] += y[bss * i + j];
1144:       }
1145:     } else {
1146:       for (i = 0; i < n; i++) {
1147:         for (j = 0; j < bss; j++) x[bs * i + idxv[j]] += y[bss * i + idxs[j]];
1148:       }
1149:     }
1150: #if !defined(PETSC_USE_COMPLEX)
1151:   } else if (addv == MAX_VALUES) {
1152:     if (!idxs) {
1153:       for (i = 0; i < n; i++) {
1154:         for (j = 0; j < bss; j++) x[bs * i + idxv[j]] = PetscMax(y[bss * i + j], x[bs * i + idxv[j]]);
1155:       }
1156:     } else {
1157:       for (i = 0; i < n; i++) {
1158:         for (j = 0; j < bss; j++) x[bs * i + idxv[j]] = PetscMax(y[bss * i + idxs[j]], x[bs * i + idxv[j]]);
1159:       }
1160:     }
1161: #endif
1162:   } else SETERRQ(PetscObjectComm((PetscObject)v), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown norm type");

1164:   PetscCall(VecRestoreArray(v, &x));
1165:   PetscCall(VecRestoreArrayRead(s, &y));
1166:   PetscFunctionReturn(PETSC_SUCCESS);
1167: }

1169: static PetscErrorCode VecApplyUnary_Private(Vec v, PetscDeviceContext dctx, const char async_name[], PetscErrorCode (*unary_op)(Vec), PetscScalar (*UnaryFunc)(PetscScalar))
1170: {
1171:   PetscFunctionBegin;
1173:   PetscCall(VecSetErrorIfLocked(v, 1));
1174:   if (dctx) {
1175:     PetscErrorCode (*unary_op_async)(Vec, PetscDeviceContext);

1177:     PetscCall(PetscObjectQueryFunction((PetscObject)v, async_name, &unary_op_async));
1178:     if (unary_op_async) {
1179:       PetscCall((*unary_op_async)(v, dctx));
1180:       PetscFunctionReturn(PETSC_SUCCESS);
1181:     }
1182:   }
1183:   if (unary_op) {
1185:     PetscCall((*unary_op)(v));
1186:   } else {
1187:     PetscInt     n;
1188:     PetscScalar *x;

1191:     PetscCall(VecGetLocalSize(v, &n));
1192:     PetscCall(VecGetArray(v, &x));
1193:     for (PetscInt i = 0; i < n; ++i) x[i] = UnaryFunc(x[i]);
1194:     PetscCall(VecRestoreArray(v, &x));
1195:   }
1196:   PetscFunctionReturn(PETSC_SUCCESS);
1197: }

1199: static PetscScalar ScalarReciprocal_Fn(PetscScalar x)
1200: {
1201:   const PetscScalar zero = 0.0;

1203:   return x == zero ? zero : ((PetscScalar)1.0) / x;
1204: }

1206: PetscErrorCode VecReciprocalAsync_Private(Vec v, PetscDeviceContext dctx)
1207: {
1208:   PetscFunctionBegin;
1209:   PetscCall(VecApplyUnary_Private(v, dctx, VecAsyncFnName(Reciprocal), v->ops->reciprocal, ScalarReciprocal_Fn));
1210:   PetscFunctionReturn(PETSC_SUCCESS);
1211: }

1213: PetscErrorCode VecReciprocal_Default(Vec v)
1214: {
1215:   PetscFunctionBegin;
1216:   PetscCall(VecApplyUnary_Private(v, NULL, NULL, NULL, ScalarReciprocal_Fn));
1217:   PetscFunctionReturn(PETSC_SUCCESS);
1218: }

1220: static PetscScalar ScalarExp_Fn(PetscScalar x)
1221: {
1222:   return PetscExpScalar(x);
1223: }

1225: PetscErrorCode VecExpAsync_Private(Vec v, PetscDeviceContext dctx)
1226: {
1227:   PetscFunctionBegin;
1229:   PetscCall(VecApplyUnary_Private(v, dctx, VecAsyncFnName(Exp), v->ops->exp, ScalarExp_Fn));
1230:   PetscFunctionReturn(PETSC_SUCCESS);
1231: }

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

1236:   Not Collective

1238:   Input Parameter:
1239: . v - The vector

1241:   Output Parameter:
1242: . v - The vector of exponents

1244:   Level: beginner

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

1248: @*/
1249: PetscErrorCode VecExp(Vec v)
1250: {
1251:   PetscFunctionBegin;
1252:   PetscCall(VecExpAsync_Private(v, NULL));
1253:   PetscFunctionReturn(PETSC_SUCCESS);
1254: }

1256: static PetscScalar ScalarLog_Fn(PetscScalar x)
1257: {
1258:   return PetscLogScalar(x);
1259: }

1261: PetscErrorCode VecLogAsync_Private(Vec v, PetscDeviceContext dctx)
1262: {
1263:   PetscFunctionBegin;
1265:   PetscCall(VecApplyUnary_Private(v, dctx, VecAsyncFnName(Log), v->ops->log, ScalarLog_Fn));
1266:   PetscFunctionReturn(PETSC_SUCCESS);
1267: }

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

1272:   Not Collective

1274:   Input Parameter:
1275: . v - The vector

1277:   Output Parameter:
1278: . v - The vector of logs

1280:   Level: beginner

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

1284: @*/
1285: PetscErrorCode VecLog(Vec v)
1286: {
1287:   PetscFunctionBegin;
1288:   PetscCall(VecLogAsync_Private(v, NULL));
1289:   PetscFunctionReturn(PETSC_SUCCESS);
1290: }

1292: static PetscScalar ScalarAbs_Fn(PetscScalar x)
1293: {
1294:   return PetscAbsScalar(x);
1295: }

1297: PetscErrorCode VecAbsAsync_Private(Vec v, PetscDeviceContext dctx)
1298: {
1299:   PetscFunctionBegin;
1301:   PetscCall(VecApplyUnary_Private(v, dctx, VecAsyncFnName(Abs), v->ops->abs, ScalarAbs_Fn));
1302:   PetscFunctionReturn(PETSC_SUCCESS);
1303: }

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

1308:   Logically Collective

1310:   Input Parameter:
1311: . v - the vector

1313:   Level: intermediate

1315: .seealso: `Vec`, `VecExp()`, `VecSqrtAbs()`, `VecReciprocal()`, `VecLog()`
1316: @*/
1317: PetscErrorCode VecAbs(Vec v)
1318: {
1319:   PetscFunctionBegin;
1320:   PetscCall(VecAbsAsync_Private(v, NULL));
1321:   PetscFunctionReturn(PETSC_SUCCESS);
1322: }

1324: static PetscScalar ScalarConjugate_Fn(PetscScalar x)
1325: {
1326:   return PetscConj(x);
1327: }

1329: PetscErrorCode VecConjugateAsync_Private(Vec v, PetscDeviceContext dctx)
1330: {
1331:   PetscFunctionBegin;
1333:   if (PetscDefined(USE_COMPLEX)) PetscCall(VecApplyUnary_Private(v, dctx, VecAsyncFnName(Conjugate), v->ops->conjugate, ScalarConjugate_Fn));
1334:   PetscFunctionReturn(PETSC_SUCCESS);
1335: }

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

1340:   Logically Collective

1342:   Input Parameter:
1343: . x - the vector

1345:   Level: intermediate

1347: .seealso: [](ch_vectors), `Vec`, `VecSet()`
1348: @*/
1349: PetscErrorCode VecConjugate(Vec x)
1350: {
1351:   PetscFunctionBegin;
1352:   PetscCall(VecConjugateAsync_Private(x, NULL));
1353:   PetscFunctionReturn(PETSC_SUCCESS);
1354: }

1356: static PetscScalar ScalarSqrtAbs_Fn(PetscScalar x)
1357: {
1358:   return PetscSqrtScalar(ScalarAbs_Fn(x));
1359: }

1361: PetscErrorCode VecSqrtAbsAsync_Private(Vec v, PetscDeviceContext dctx)
1362: {
1363:   PetscFunctionBegin;
1365:   PetscCall(VecApplyUnary_Private(v, dctx, VecAsyncFnName(SqrtAbs), v->ops->sqrt, ScalarSqrtAbs_Fn));
1366:   PetscFunctionReturn(PETSC_SUCCESS);
1367: }

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

1372:   Not Collective

1374:   Input Parameter:
1375: . v - The vector

1377:   Level: beginner

1379:   Note:
1380:   The actual function is sqrt(|x_i|)

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

1384: @*/
1385: PetscErrorCode VecSqrtAbs(Vec v)
1386: {
1387:   PetscFunctionBegin;
1388:   PetscCall(VecSqrtAbsAsync_Private(v, NULL));
1389:   PetscFunctionReturn(PETSC_SUCCESS);
1390: }

1392: static PetscScalar ScalarImaginaryPart_Fn(PetscScalar x)
1393: {
1394:   const PetscReal imag = PetscImaginaryPart(x);

1396: #if PetscDefined(USE_COMPLEX)
1397:   return PetscCMPLX(imag, 0.0);
1398: #else
1399:   return imag;
1400: #endif
1401: }

1403: /*@
1404:   VecImaginaryPart - Replaces a complex vector with its imginary part

1406:   Collective

1408:   Input Parameter:
1409: . v - the vector

1411:   Level: beginner

1413: .seealso: `Vec`, `VecNorm()`, `VecRealPart()`
1414: @*/
1415: PetscErrorCode VecImaginaryPart(Vec v)
1416: {
1417:   PetscFunctionBegin;
1419:   PetscCall(VecApplyUnary_Private(v, NULL, NULL, NULL, ScalarImaginaryPart_Fn));
1420:   PetscFunctionReturn(PETSC_SUCCESS);
1421: }

1423: static PetscScalar ScalarRealPart_Fn(PetscScalar x)
1424: {
1425:   const PetscReal real = PetscRealPart(x);

1427: #if PetscDefined(USE_COMPLEX)
1428:   return PetscCMPLX(real, 0.0);
1429: #else
1430:   return real;
1431: #endif
1432: }

1434: /*@
1435:   VecRealPart - Replaces a complex vector with its real part

1437:   Collective

1439:   Input Parameter:
1440: . v - the vector

1442:   Level: beginner

1444: .seealso: `Vec`, `VecNorm()`, `VecImaginaryPart()`
1445: @*/
1446: PetscErrorCode VecRealPart(Vec v)
1447: {
1448:   PetscFunctionBegin;
1450:   PetscCall(VecApplyUnary_Private(v, NULL, NULL, NULL, ScalarRealPart_Fn));
1451:   PetscFunctionReturn(PETSC_SUCCESS);
1452: }

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

1457:   Collective

1459:   Input Parameters:
1460: + s - first vector
1461: - t - second vector

1463:   Output Parameters:
1464: + dp - s'conj(t)
1465: - nm - t'conj(t)

1467:   Level: advanced

1469:   Note:
1470:   conj(x) is the complex conjugate of x when x is complex

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

1474: @*/
1475: PetscErrorCode VecDotNorm2(Vec s, Vec t, PetscScalar *dp, PetscReal *nm)
1476: {
1477:   PetscScalar work[] = {0.0, 0.0};

1479:   PetscFunctionBegin;
1482:   PetscAssertPointer(dp, 3);
1483:   PetscAssertPointer(nm, 4);
1486:   PetscCheckSameTypeAndComm(s, 1, t, 2);
1487:   PetscCheck(s->map->N == t->map->N, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Incompatible vector global lengths");
1488:   PetscCheck(s->map->n == t->map->n, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Incompatible vector local lengths");

1490:   PetscCall(PetscLogEventBegin(VEC_DotNorm2, s, t, 0, 0));
1491:   if (s->ops->dotnorm2) {
1492:     PetscUseTypeMethod(s, dotnorm2, t, work, work + 1);
1493:   } else {
1494:     const PetscScalar *sx, *tx;
1495:     PetscInt           n;

1497:     PetscCall(VecGetLocalSize(s, &n));
1498:     PetscCall(VecGetArrayRead(s, &sx));
1499:     PetscCall(VecGetArrayRead(t, &tx));
1500:     for (PetscInt i = 0; i < n; ++i) {
1501:       const PetscScalar txconj = PetscConj(tx[i]);

1503:       work[0] += sx[i] * txconj;
1504:       work[1] += tx[i] * txconj;
1505:     }
1506:     PetscCall(VecRestoreArrayRead(t, &tx));
1507:     PetscCall(VecRestoreArrayRead(s, &sx));
1508:     PetscCall(MPIU_Allreduce(MPI_IN_PLACE, work, 2, MPIU_SCALAR, MPIU_SUM, PetscObjectComm((PetscObject)s)));
1509:     PetscCall(PetscLogFlops(4.0 * n));
1510:   }
1511:   PetscCall(PetscLogEventEnd(VEC_DotNorm2, s, t, 0, 0));
1512:   *dp = work[0];
1513:   *nm = PetscRealPart(work[1]);
1514:   PetscFunctionReturn(PETSC_SUCCESS);
1515: }

1517: /*@
1518:   VecSum - Computes the sum of all the components of a vector.

1520:   Collective

1522:   Input Parameter:
1523: . v - the vector

1525:   Output Parameter:
1526: . sum - the result

1528:   Level: beginner

1530: .seealso: `Vec`, `VecMean()`, `VecNorm()`
1531: @*/
1532: PetscErrorCode VecSum(Vec v, PetscScalar *sum)
1533: {
1534:   PetscScalar tmp = 0.0;

1536:   PetscFunctionBegin;
1538:   PetscAssertPointer(sum, 2);
1539:   if (v->ops->sum) {
1540:     PetscUseTypeMethod(v, sum, &tmp);
1541:   } else {
1542:     const PetscScalar *x;
1543:     PetscInt           n;

1545:     PetscCall(VecGetLocalSize(v, &n));
1546:     PetscCall(VecGetArrayRead(v, &x));
1547:     for (PetscInt i = 0; i < n; ++i) tmp += x[i];
1548:     PetscCall(VecRestoreArrayRead(v, &x));
1549:   }
1550:   PetscCall(MPIU_Allreduce(MPI_IN_PLACE, &tmp, 1, MPIU_SCALAR, MPIU_SUM, PetscObjectComm((PetscObject)v)));
1551:   *sum = tmp;
1552:   PetscFunctionReturn(PETSC_SUCCESS);
1553: }

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

1558:   Collective

1560:   Input Parameter:
1561: . v - the vector

1563:   Output Parameter:
1564: . mean - the result

1566:   Level: beginner

1568: .seealso: `Vec`, `VecSum()`, `VecNorm()`
1569: @*/
1570: PetscErrorCode VecMean(Vec v, PetscScalar *mean)
1571: {
1572:   PetscInt n;

1574:   PetscFunctionBegin;
1576:   PetscAssertPointer(mean, 2);
1577:   PetscCall(VecGetSize(v, &n));
1578:   PetscCall(VecSum(v, mean));
1579:   *mean /= n;
1580:   PetscFunctionReturn(PETSC_SUCCESS);
1581: }

1583: PetscErrorCode VecShiftAsync_Private(Vec v, PetscScalar shift, PetscDeviceContext dctx)
1584: {
1585:   PetscErrorCode (*shift_async)(Vec, PetscScalar, PetscDeviceContext) = NULL;

1587:   PetscFunctionBegin;
1590:   PetscCall(VecSetErrorIfLocked(v, 1));
1591:   if (shift == (PetscScalar)0.0) PetscFunctionReturn(PETSC_SUCCESS);

1593:   if (dctx) {
1594:     PetscErrorCode (*shift_async)(Vec, PetscScalar, PetscDeviceContext);

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

1606:     PetscCall(VecGetLocalSize(v, &n));
1607:     PetscCall(VecGetArray(v, &x));
1608:     for (PetscInt i = 0; i < n; ++i) x[i] += shift;
1609:     PetscCall(VecRestoreArray(v, &x));
1610:     PetscCall(PetscLogFlops(n));
1611:   }
1612:   PetscFunctionReturn(PETSC_SUCCESS);
1613: }

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

1619:   Logically Collective

1621:   Input Parameters:
1622: + v     - the vector
1623: - shift - the shift

1625:   Level: intermediate

1627: .seealso: `Vec`
1628: @*/
1629: PetscErrorCode VecShift(Vec v, PetscScalar shift)
1630: {
1631:   PetscFunctionBegin;
1632:   PetscCall(VecShiftAsync_Private(v, shift, NULL));
1633:   PetscFunctionReturn(PETSC_SUCCESS);
1634: }

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

1639:   Input Parameters:
1640: + x   - The vector
1641: . row - The ordering
1642: - inv - The flag for inverting the permutation

1644:   Level: beginner

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

1649: .seealso: `Vec`, `MatPermute()`
1650: @*/
1651: PetscErrorCode VecPermute(Vec x, IS row, PetscBool inv)
1652: {
1653:   const PetscScalar *array;
1654:   PetscScalar       *newArray;
1655:   const PetscInt    *idx;
1656:   PetscInt           i, rstart, rend;

1658:   PetscFunctionBegin;
1661:   PetscCall(VecSetErrorIfLocked(x, 1));
1662:   PetscCall(VecGetOwnershipRange(x, &rstart, &rend));
1663:   PetscCall(ISGetIndices(row, &idx));
1664:   PetscCall(VecGetArrayRead(x, &array));
1665:   PetscCall(PetscMalloc1(x->map->n, &newArray));
1666:   if (PetscDefined(USE_DEBUG)) {
1667:     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]);
1668:   }
1669:   if (!inv) {
1670:     for (i = 0; i < x->map->n; i++) newArray[i] = array[idx[i] - rstart];
1671:   } else {
1672:     for (i = 0; i < x->map->n; i++) newArray[idx[i] - rstart] = array[i];
1673:   }
1674:   PetscCall(VecRestoreArrayRead(x, &array));
1675:   PetscCall(ISRestoreIndices(row, &idx));
1676:   PetscCall(VecReplaceArray(x, newArray));
1677:   PetscFunctionReturn(PETSC_SUCCESS);
1678: }

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

1685:   Collective

1687:   Input Parameters:
1688: + vec1 - the first vector
1689: - vec2 - the second vector

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

1694:   Level: intermediate

1696: .seealso: `Vec`
1697: @*/
1698: PetscErrorCode VecEqual(Vec vec1, Vec vec2, PetscBool *flg)
1699: {
1700:   const PetscScalar *v1, *v2;
1701:   PetscInt           n1, n2, N1, N2;
1702:   PetscBool          flg1;

1704:   PetscFunctionBegin;
1707:   PetscAssertPointer(flg, 3);
1708:   if (vec1 == vec2) *flg = PETSC_TRUE;
1709:   else {
1710:     PetscCall(VecGetSize(vec1, &N1));
1711:     PetscCall(VecGetSize(vec2, &N2));
1712:     if (N1 != N2) flg1 = PETSC_FALSE;
1713:     else {
1714:       PetscCall(VecGetLocalSize(vec1, &n1));
1715:       PetscCall(VecGetLocalSize(vec2, &n2));
1716:       if (n1 != n2) flg1 = PETSC_FALSE;
1717:       else {
1718:         PetscCall(VecGetArrayRead(vec1, &v1));
1719:         PetscCall(VecGetArrayRead(vec2, &v2));
1720:         PetscCall(PetscArraycmp(v1, v2, n1, &flg1));
1721:         PetscCall(VecRestoreArrayRead(vec1, &v1));
1722:         PetscCall(VecRestoreArrayRead(vec2, &v2));
1723:       }
1724:     }
1725:     /* combine results from all processors */
1726:     PetscCall(MPIU_Allreduce(&flg1, flg, 1, MPIU_BOOL, MPI_MIN, PetscObjectComm((PetscObject)vec1)));
1727:   }
1728:   PetscFunctionReturn(PETSC_SUCCESS);
1729: }

1731: /*@
1732:   VecUniqueEntries - Compute the number of unique entries, and those entries

1734:   Collective

1736:   Input Parameter:
1737: . vec - the vector

1739:   Output Parameters:
1740: + n - The number of unique entries
1741: - e - The entries

1743:   Level: intermediate

1745: .seealso: `Vec`
1746: @*/
1747: PetscErrorCode VecUniqueEntries(Vec vec, PetscInt *n, PetscScalar **e)
1748: {
1749:   const PetscScalar *v;
1750:   PetscScalar       *tmp, *vals;
1751:   PetscMPIInt       *N, *displs, l;
1752:   PetscInt           ng, m, i, j, p;
1753:   PetscMPIInt        size;

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