Actual source code: vinv.c


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

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

 11:    Logically Collective

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

 18:    Level: advanced

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

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

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

 33:   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\n  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()`, `VecStrideScale()`
 65: @*/
 66: PetscErrorCode VecStrideScale(Vec v, PetscInt start, PetscScalar scale)
 67: {
 68:   PetscInt     i, n, bs;
 69:   PetscScalar *x;

 71:   PetscFunctionBegin;
 73:   PetscCall(VecGetLocalSize(v, &n));
 74:   PetscCall(VecGetBlockSize(v, &bs));
 75:   PetscCheck(start >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Negative start %" PetscInt_FMT, start);
 76:   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);
 77:   PetscCall(VecGetArray(v, &x));
 78:   for (i = start; i < n; i += bs) x[i] *= scale;
 79:   PetscCall(VecRestoreArray(v, &x));
 80:   PetscFunctionReturn(PETSC_SUCCESS);
 81: }

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

 87:    Collective

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

 94:    Output Parameter:
 95: .  norm - the norm

 97:    Level: advanced

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

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

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

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

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

119:   PetscFunctionBegin;
123:   PetscCall(VecGetLocalSize(v, &n));
124:   PetscCall(VecGetBlockSize(v, &bs));
125:   PetscCheck(start >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Negative start %" PetscInt_FMT, start);
126:   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);
127:   PetscCall(VecGetArrayRead(v, &x));
128:   if (ntype == NORM_2) {
129:     PetscScalar sum = 0.0;
130:     for (i = start; i < n; i += bs) sum += x[i] * (PetscConj(x[i]));
131:     tnorm = PetscRealPart(sum);
132:     PetscCall(MPIU_Allreduce(&tnorm, nrm, 1, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)v)));
133:     *nrm = PetscSqrtReal(*nrm);
134:   } else if (ntype == NORM_1) {
135:     tnorm = 0.0;
136:     for (i = start; i < n; i += bs) tnorm += PetscAbsScalar(x[i]);
137:     PetscCall(MPIU_Allreduce(&tnorm, nrm, 1, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)v)));
138:   } else if (ntype == NORM_INFINITY) {
139:     tnorm = 0.0;
140:     for (i = start; i < n; i += bs) {
141:       if (PetscAbsScalar(x[i]) > tnorm) tnorm = PetscAbsScalar(x[i]);
142:     }
143:     PetscCall(MPIU_Allreduce(&tnorm, nrm, 1, MPIU_REAL, MPIU_MAX, PetscObjectComm((PetscObject)v)));
144:   } else SETERRQ(PetscObjectComm((PetscObject)v), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown norm type");
145:   PetscCall(VecRestoreArrayRead(v, &x));
146:   PetscFunctionReturn(PETSC_SUCCESS);
147: }

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

153:    Collective

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

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

163:    Level: advanced

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

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

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

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

184:   PetscFunctionBegin;
187:   PetscCall(VecGetLocalSize(v, &n));
188:   PetscCall(VecGetBlockSize(v, &bs));
189:   PetscCheck(start >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Negative start %" PetscInt_FMT, start);
190:   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);
191:   PetscCall(VecGetArrayRead(v, &x));
192:   for (i = start; i < n; i += bs) {
193:     if (PetscRealPart(x[i]) > max) {
194:       max = PetscRealPart(x[i]);
195:       id  = i;
196:     }
197:   }
198:   PetscCall(VecRestoreArrayRead(v, &x));
199: #if defined(PETSC_HAVE_MPIUNI)
200:   *nrm = max;
201:   if (idex) *idex = id;
202: #else
203:   if (!idex) {
204:     PetscCall(MPIU_Allreduce(&max, nrm, 1, MPIU_REAL, MPIU_MAX, PetscObjectComm((PetscObject)v)));
205:   } else {
206:     struct {
207:       PetscReal v;
208:       PetscInt  i;
209:     } in, out;
210:     PetscInt rstart;

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

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

227:    Collective

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

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

237:    Level: advanced

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

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

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

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

258:   PetscFunctionBegin;
261:   PetscCall(VecGetLocalSize(v, &n));
262:   PetscCall(VecGetBlockSize(v, &bs));
263:   PetscCheck(start >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Negative start %" PetscInt_FMT, start);
264:   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);
265:   PetscCall(VecGetArrayRead(v, &x));
266:   for (i = start; i < n; i += bs) {
267:     if (PetscRealPart(x[i]) < min) {
268:       min = PetscRealPart(x[i]);
269:       id  = i;
270:     }
271:   }
272:   PetscCall(VecRestoreArrayRead(v, &x));
273: #if defined(PETSC_HAVE_MPIUNI)
274:   *nrm = min;
275:   if (idex) *idex = id;
276: #else
277:   if (!idex) {
278:     PetscCall(MPIU_Allreduce(&min, nrm, 1, MPIU_REAL, MPIU_MIN, PetscObjectComm((PetscObject)v)));
279:   } else {
280:     struct {
281:       PetscReal v;
282:       PetscInt  i;
283:     } in, out;
284:     PetscInt rstart;

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

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

301:    Collective

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

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

310:    Level: advanced

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

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

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

327:   PetscFunctionBegin;
330:   PetscCall(VecGetLocalSize(v, &n));
331:   PetscCall(VecGetBlockSize(v, &bs));
332:   PetscCheck(start >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Negative start %" PetscInt_FMT, start);
333:   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);
334:   PetscCall(VecGetArrayRead(v, &x));
335:   for (i = start; i < n; i += bs) local_sum += x[i];
336:   PetscCall(MPIU_Allreduce(&local_sum, sum, 1, MPIU_SCALAR, MPIU_SUM, PetscObjectComm((PetscObject)v)));
337:   PetscCall(VecRestoreArrayRead(v, &x));
338:   PetscFunctionReturn(PETSC_SUCCESS);
339: }

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

345:    Logically Collective

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

351:    Level: advanced

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

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

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

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

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

384:    Collective

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

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

393:    Level: advanced

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

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

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

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

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

415:   PetscFunctionBegin;
419:   PetscCall(VecGetLocalSize(v, &n));
420:   PetscCall(VecGetArrayRead(v, &x));
421:   PetscCall(PetscObjectGetComm((PetscObject)v, &comm));

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

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

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

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

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

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

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

468:    Collective

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

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

478:    Level: advanced

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

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

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

495:   PetscFunctionBegin;
498:   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");
499:   PetscCall(VecGetLocalSize(v, &n));
500:   PetscCall(VecGetArrayRead(v, &x));
501:   PetscCall(PetscObjectGetComm((PetscObject)v, &comm));

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

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

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

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

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

527:    Collective

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

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

537:    Level: advanced

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

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

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

554:   PetscFunctionBegin;
557:   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");
558:   PetscCall(VecGetLocalSize(v, &n));
559:   PetscCall(VecGetArrayRead(v, &x));
560:   PetscCall(PetscObjectGetComm((PetscObject)v, &comm));

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

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

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

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

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

585:    Collective

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

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

593:    Level: advanced

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

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

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

611:   PetscFunctionBegin;
614:   PetscCall(VecGetLocalSize(v, &n));
615:   PetscCall(VecGetArrayRead(v, &x));
616:   PetscCall(PetscObjectGetComm((PetscObject)v, &comm));

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

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

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

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

636:    Collective

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

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

645:    Level: advanced

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

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

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

658:    Not optimized; could be easily

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

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

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

692:   n = n / bs;

694:   jj = 0;
695:   if (addv == INSERT_VALUES) {
696:     for (j = 0; j < nv; j++) {
697:       for (k = 0; k < bss[j]; k++) {
698:         for (i = 0; i < n; i++) y[j][i * bss[j] + k] = x[bs * i + jj + k];
699:       }
700:       jj += bss[j];
701:     }
702:   } else if (addv == ADD_VALUES) {
703:     for (j = 0; j < nv; j++) {
704:       for (k = 0; k < bss[j]; k++) {
705:         for (i = 0; i < n; i++) y[j][i * bss[j] + k] += x[bs * i + jj + k];
706:       }
707:       jj += bss[j];
708:     }
709: #if !defined(PETSC_USE_COMPLEX)
710:   } else if (addv == MAX_VALUES) {
711:     for (j = 0; j < nv; j++) {
712:       for (k = 0; k < bss[j]; k++) {
713:         for (i = 0; i < n; i++) y[j][i * bss[j] + k] = PetscMax(y[j][i * bss[j] + k], x[bs * i + jj + k]);
714:       }
715:       jj += bss[j];
716:     }
717: #endif
718:   } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown insert type");

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

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

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

731:    Collective

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

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

740:    Level: advanced

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

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

749:    Not optimized; could be easily

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

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

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

783:   n = n / bs;

785:   jj = 0;
786:   if (addv == INSERT_VALUES) {
787:     for (j = 0; j < nv; j++) {
788:       for (k = 0; k < bss[j]; k++) {
789:         for (i = 0; i < n; i++) x[bs * i + jj + k] = y[j][i * bss[j] + k];
790:       }
791:       jj += bss[j];
792:     }
793:   } else if (addv == ADD_VALUES) {
794:     for (j = 0; j < nv; j++) {
795:       for (k = 0; k < bss[j]; k++) {
796:         for (i = 0; i < n; i++) x[bs * i + jj + k] += y[j][i * bss[j] + k];
797:       }
798:       jj += bss[j];
799:     }
800: #if !defined(PETSC_USE_COMPLEX)
801:   } else if (addv == MAX_VALUES) {
802:     for (j = 0; j < nv; j++) {
803:       for (k = 0; k < bss[j]; k++) {
804:         for (i = 0; i < n; i++) x[bs * i + jj + k] = PetscMax(x[bs * i + jj + k], y[j][i * bss[j] + k]);
805:       }
806:       jj += bss[j];
807:     }
808: #endif
809:   } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown insert type");

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

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

821:    Collective

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

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

831:    Level: advanced

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

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

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

843:    Not optimized; could be easily

845: .seealso: `Vec`, `VecStrideNorm()`, `VecStrideScatter()`, `VecStrideMin()`, `VecStrideMax()`, `VecStrideGatherAll()`,
846:           `VecStrideScatterAll()`
847: @*/
848: PetscErrorCode VecStrideGather(Vec v, PetscInt start, Vec s, InsertMode addv)
849: {
850:   PetscFunctionBegin;
853:   PetscCheck(start >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Negative start %" PetscInt_FMT, start);
854:   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,
855:              v->map->bs);
856:   PetscUseTypeMethod(v, stridegather, start, s, addv);
857:   PetscFunctionReturn(PETSC_SUCCESS);
858: }

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

863:    Collective

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

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

873:    Level: advanced

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

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

882:    Not optimized; could be easily

884: .seealso: `Vec`, `VecStrideNorm()`, `VecStrideGather()`, `VecStrideMin()`, `VecStrideMax()`, `VecStrideGatherAll()`,
885:           `VecStrideScatterAll()`, `VecStrideSubSetScatter()`, `VecStrideSubSetGather()`
886: @*/
887: PetscErrorCode VecStrideScatter(Vec s, PetscInt start, Vec v, InsertMode addv)
888: {
889:   PetscFunctionBegin;
892:   PetscCheck(start >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Negative start %" PetscInt_FMT, start);
893:   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,
894:              v->map->bs);
895:   PetscCall((*v->ops->stridescatter)(s, start, v, addv));
896:   PetscFunctionReturn(PETSC_SUCCESS);
897: }

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

903:    Collective

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

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

915:    Level: advanced

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

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

923:    Not optimized; could be easily

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

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

941:    Collective

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

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

953:    Level: advanced

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

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

961:    Not optimized; could be easily

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

1170: static PetscErrorCode VecApplyUnary_Private(Vec v, PetscErrorCode (*unary_op)(Vec), PetscScalar (*UnaryFunc)(PetscScalar))
1171: {
1172:   PetscFunctionBegin;
1174:   PetscCall(VecSetErrorIfLocked(v, 1));
1175:   if (unary_op) {
1177:     PetscCall((*unary_op)(v));
1178:   } else {
1179:     PetscInt     n;
1180:     PetscScalar *x;

1183:     PetscCall(VecGetLocalSize(v, &n));
1184:     PetscCall(VecGetArray(v, &x));
1185:     for (PetscInt i = 0; i < n; ++i) x[i] = UnaryFunc(x[i]);
1186:     PetscCall(VecRestoreArray(v, &x));
1187:   }
1188:   PetscFunctionReturn(PETSC_SUCCESS);
1189: }

1191: static PetscScalar ScalarReciprocal_Fn(PetscScalar x)
1192: {
1193:   const PetscScalar zero = 0.0;

1195:   return x == zero ? zero : ((PetscScalar)1.0) / x;
1196: }

1198: PetscErrorCode VecReciprocal_Default(Vec v)
1199: {
1200:   PetscFunctionBegin;
1201:   PetscCall(VecApplyUnary_Private(v, NULL, ScalarReciprocal_Fn));
1202:   PetscFunctionReturn(PETSC_SUCCESS);
1203: }

1205: static PetscScalar ScalarExp_Fn(PetscScalar x)
1206: {
1207:   return PetscExpScalar(x);
1208: }

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

1213:   Not Collective

1215:   Input Parameter:
1216: . v - The vector

1218:   Output Parameter:
1219: . v - The vector of exponents

1221:   Level: beginner

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

1225: @*/
1226: PetscErrorCode VecExp(Vec v)
1227: {
1228:   PetscFunctionBegin;
1230:   PetscCall(VecApplyUnary_Private(v, v->ops->exp, ScalarExp_Fn));
1231:   PetscFunctionReturn(PETSC_SUCCESS);
1232: }

1234: static PetscScalar ScalarLog_Fn(PetscScalar x)
1235: {
1236:   return PetscLogScalar(x);
1237: }

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

1242:   Not Collective

1244:   Input Parameter:
1245: . v - The vector

1247:   Output Parameter:
1248: . v - The vector of logs

1250:   Level: beginner

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

1254: @*/
1255: PetscErrorCode VecLog(Vec v)
1256: {
1257:   PetscFunctionBegin;
1259:   PetscCall(VecApplyUnary_Private(v, v->ops->log, ScalarLog_Fn));
1260:   PetscFunctionReturn(PETSC_SUCCESS);
1261: }

1263: static PetscScalar ScalarAbs_Fn(PetscScalar x)
1264: {
1265:   return PetscAbsScalar(x);
1266: }

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

1271:    Logically Collective

1273:    Input Parameter:
1274: .  v - the vector

1276:    Level: intermediate

1278: @*/
1279: PetscErrorCode VecAbs(Vec v)
1280: {
1281:   PetscFunctionBegin;
1283:   PetscCall(VecApplyUnary_Private(v, v->ops->abs, ScalarAbs_Fn));
1284:   PetscFunctionReturn(PETSC_SUCCESS);
1285: }

1287: static PetscScalar ScalarSqrtAbs_Fn(PetscScalar x)
1288: {
1289:   return PetscSqrtScalar(ScalarAbs_Fn(x));
1290: }

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

1295:   Not Collective

1297:   Input Parameter:
1298: . v - The vector

1300:   Level: beginner

1302:   Note:
1303:   The actual function is sqrt(|x_i|)

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

1307: @*/
1308: PetscErrorCode VecSqrtAbs(Vec v)
1309: {
1310:   PetscFunctionBegin;
1312:   PetscCall(VecApplyUnary_Private(v, v->ops->sqrt, ScalarSqrtAbs_Fn));
1313:   PetscFunctionReturn(PETSC_SUCCESS);
1314: }

1316: static PetscScalar ScalarImaginaryPart_Fn(PetscScalar x)
1317: {
1318:   const PetscReal imag = PetscImaginaryPart(x);

1320: #if PetscDefined(USE_COMPLEX)
1321:   return PetscCMPLX(imag, 0.0);
1322: #else
1323:   return imag;
1324: #endif
1325: }

1327: /*@
1328:    VecImaginaryPart - Replaces a complex vector with its imginary part

1330:    Collective

1332:    Input Parameter:
1333: .  v - the vector

1335:    Level: beginner

1337: .seealso: `Vec`, `VecNorm()`, `VecRealPart()`
1338: @*/
1339: PetscErrorCode VecImaginaryPart(Vec v)
1340: {
1341:   PetscFunctionBegin;
1343:   PetscCall(VecApplyUnary_Private(v, NULL, ScalarImaginaryPart_Fn));
1344:   PetscFunctionReturn(PETSC_SUCCESS);
1345: }

1347: static PetscScalar ScalarRealPart_Fn(PetscScalar x)
1348: {
1349:   const PetscReal real = PetscRealPart(x);

1351: #if PetscDefined(USE_COMPLEX)
1352:   return PetscCMPLX(real, 0.0);
1353: #else
1354:   return real;
1355: #endif
1356: }

1358: /*@
1359:    VecRealPart - Replaces a complex vector with its real part

1361:    Collective

1363:    Input Parameter:
1364: .  v - the vector

1366:    Level: beginner

1368: .seealso: `Vec`, `VecNorm()`, `VecImaginaryPart()`
1369: @*/
1370: PetscErrorCode VecRealPart(Vec v)
1371: {
1372:   PetscFunctionBegin;
1374:   PetscCall(VecApplyUnary_Private(v, NULL, ScalarRealPart_Fn));
1375:   PetscFunctionReturn(PETSC_SUCCESS);
1376: }

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

1381:   Collective

1383:   Input Parameters:
1384: + s - first vector
1385: - t - second vector

1387:   Output Parameters:
1388: + dp - s'conj(t)
1389: - nm - t'conj(t)

1391:   Level: advanced

1393:   Note:
1394:     conj(x) is the complex conjugate of x when x is complex

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

1398: @*/
1399: PetscErrorCode VecDotNorm2(Vec s, Vec t, PetscScalar *dp, PetscReal *nm)
1400: {
1401:   PetscScalar work[] = {0.0, 0.0};

1403:   PetscFunctionBegin;
1410:   PetscCheckSameTypeAndComm(s, 1, t, 2);
1411:   PetscCheck(s->map->N == t->map->N, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Incompatible vector global lengths");
1412:   PetscCheck(s->map->n == t->map->n, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Incompatible vector local lengths");

1414:   PetscCall(PetscLogEventBegin(VEC_DotNorm2, s, t, 0, 0));
1415:   if (s->ops->dotnorm2) {
1416:     PetscUseTypeMethod(s, dotnorm2, t, work, work + 1);
1417:   } else {
1418:     const PetscScalar *sx, *tx;
1419:     PetscInt           n;

1421:     PetscCall(VecGetLocalSize(s, &n));
1422:     PetscCall(VecGetArrayRead(s, &sx));
1423:     PetscCall(VecGetArrayRead(t, &tx));
1424:     for (PetscInt i = 0; i < n; ++i) {
1425:       const PetscScalar txconj = PetscConj(tx[i]);

1427:       work[0] += sx[i] * txconj;
1428:       work[1] += tx[i] * txconj;
1429:     }
1430:     PetscCall(VecRestoreArrayRead(t, &tx));
1431:     PetscCall(VecRestoreArrayRead(s, &sx));
1432:     PetscCall(MPIU_Allreduce(MPI_IN_PLACE, work, 2, MPIU_SCALAR, MPIU_SUM, PetscObjectComm((PetscObject)s)));
1433:     PetscCall(PetscLogFlops(4.0 * n));
1434:   }
1435:   PetscCall(PetscLogEventEnd(VEC_DotNorm2, s, t, 0, 0));
1436:   *dp = work[0];
1437:   *nm = PetscRealPart(work[1]);
1438:   PetscFunctionReturn(PETSC_SUCCESS);
1439: }

1441: /*@
1442:    VecSum - Computes the sum of all the components of a vector.

1444:    Collective

1446:    Input Parameter:
1447: .  v - the vector

1449:    Output Parameter:
1450: .  sum - the result

1452:    Level: beginner

1454: .seealso: `Vec`, `VecMean()`, `VecNorm()`
1455: @*/
1456: PetscErrorCode VecSum(Vec v, PetscScalar *sum)
1457: {
1458:   PetscScalar tmp = 0.0;

1460:   PetscFunctionBegin;
1463:   if (v->ops->sum) {
1464:     PetscUseTypeMethod(v, sum, &tmp);
1465:   } else {
1466:     const PetscScalar *x;
1467:     PetscInt           n;

1469:     PetscCall(VecGetLocalSize(v, &n));
1470:     PetscCall(VecGetArrayRead(v, &x));
1471:     for (PetscInt i = 0; i < n; ++i) tmp += x[i];
1472:     PetscCall(VecRestoreArrayRead(v, &x));
1473:   }
1474:   PetscCall(MPIU_Allreduce(MPI_IN_PLACE, &tmp, 1, MPIU_SCALAR, MPIU_SUM, PetscObjectComm((PetscObject)v)));
1475:   *sum = tmp;
1476:   PetscFunctionReturn(PETSC_SUCCESS);
1477: }

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

1482:    Collective

1484:    Input Parameter:
1485: .  v - the vector

1487:    Output Parameter:
1488: .  mean - the result

1490:    Level: beginner

1492: .seealso: `Vec`, `VecSum()`, `VecNorm()`
1493: @*/
1494: PetscErrorCode VecMean(Vec v, PetscScalar *mean)
1495: {
1496:   PetscInt n;

1498:   PetscFunctionBegin;
1501:   PetscCall(VecGetSize(v, &n));
1502:   PetscCall(VecSum(v, mean));
1503:   *mean /= n;
1504:   PetscFunctionReturn(PETSC_SUCCESS);
1505: }

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

1511:    Logically Collective

1513:    Input Parameters:
1514: +  v - the vector
1515: -  shift - the shift

1517:    Level: intermediate

1519: .seealso: `Vec`
1520: @*/
1521: PetscErrorCode VecShift(Vec v, PetscScalar shift)
1522: {
1523:   PetscFunctionBegin;
1526:   PetscCall(VecSetErrorIfLocked(v, 1));
1527:   if (shift == (PetscScalar)0.0) PetscFunctionReturn(PETSC_SUCCESS);

1529:   if (v->ops->shift) {
1530:     PetscUseTypeMethod(v, shift, shift);
1531:   } else {
1532:     PetscInt     n;
1533:     PetscScalar *x;

1535:     PetscCall(VecGetLocalSize(v, &n));
1536:     PetscCall(VecGetArray(v, &x));
1537:     for (PetscInt i = 0; i < n; ++i) x[i] += shift;
1538:     PetscCall(VecRestoreArray(v, &x));
1539:   }
1540:   PetscFunctionReturn(PETSC_SUCCESS);
1541: }

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

1546:   Input Parameters:
1547: + vec   - The vector
1548: . order - The ordering
1549: - inv   - The flag for inverting the permutation

1551:   Level: beginner

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

1556: .seealso: `Vec`, `MatPermute()`
1557: @*/
1558: PetscErrorCode VecPermute(Vec x, IS row, PetscBool inv)
1559: {
1560:   const PetscScalar *array;
1561:   PetscScalar       *newArray;
1562:   const PetscInt    *idx;
1563:   PetscInt           i, rstart, rend;

1565:   PetscFunctionBegin;
1568:   PetscCall(VecSetErrorIfLocked(x, 1));
1569:   PetscCall(VecGetOwnershipRange(x, &rstart, &rend));
1570:   PetscCall(ISGetIndices(row, &idx));
1571:   PetscCall(VecGetArrayRead(x, &array));
1572:   PetscCall(PetscMalloc1(x->map->n, &newArray));
1573:   if (PetscDefined(USE_DEBUG)) {
1574:     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]);
1575:   }
1576:   if (!inv) {
1577:     for (i = 0; i < x->map->n; i++) newArray[i] = array[idx[i] - rstart];
1578:   } else {
1579:     for (i = 0; i < x->map->n; i++) newArray[idx[i] - rstart] = array[i];
1580:   }
1581:   PetscCall(VecRestoreArrayRead(x, &array));
1582:   PetscCall(ISRestoreIndices(row, &idx));
1583:   PetscCall(VecReplaceArray(x, newArray));
1584:   PetscFunctionReturn(PETSC_SUCCESS);
1585: }

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

1592:    Collective

1594:    Input Parameters:
1595: +  vec1 - the first vector
1596: -  vec2 - the second vector

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

1601:    Level: intermediate

1603: .seealso: `Vec`
1604: @*/
1605: PetscErrorCode VecEqual(Vec vec1, Vec vec2, PetscBool *flg)
1606: {
1607:   const PetscScalar *v1, *v2;
1608:   PetscInt           n1, n2, N1, N2;
1609:   PetscBool          flg1;

1611:   PetscFunctionBegin;
1615:   if (vec1 == vec2) *flg = PETSC_TRUE;
1616:   else {
1617:     PetscCall(VecGetSize(vec1, &N1));
1618:     PetscCall(VecGetSize(vec2, &N2));
1619:     if (N1 != N2) flg1 = PETSC_FALSE;
1620:     else {
1621:       PetscCall(VecGetLocalSize(vec1, &n1));
1622:       PetscCall(VecGetLocalSize(vec2, &n2));
1623:       if (n1 != n2) flg1 = PETSC_FALSE;
1624:       else {
1625:         PetscCall(VecGetArrayRead(vec1, &v1));
1626:         PetscCall(VecGetArrayRead(vec2, &v2));
1627:         PetscCall(PetscArraycmp(v1, v2, n1, &flg1));
1628:         PetscCall(VecRestoreArrayRead(vec1, &v1));
1629:         PetscCall(VecRestoreArrayRead(vec2, &v2));
1630:       }
1631:     }
1632:     /* combine results from all processors */
1633:     PetscCall(MPIU_Allreduce(&flg1, flg, 1, MPIU_BOOL, MPI_MIN, PetscObjectComm((PetscObject)vec1)));
1634:   }
1635:   PetscFunctionReturn(PETSC_SUCCESS);
1636: }

1638: /*@
1639:    VecUniqueEntries - Compute the number of unique entries, and those entries

1641:    Collective

1643:    Input Parameter:
1644: .  vec - the vector

1646:    Output Parameters:
1647: +  n - The number of unique entries
1648: -  e - The entries

1650:    Level: intermediate

1652: .seealso: `Vec`
1653: @*/
1654: PetscErrorCode VecUniqueEntries(Vec vec, PetscInt *n, PetscScalar **e)
1655: {
1656:   const PetscScalar *v;
1657:   PetscScalar       *tmp, *vals;
1658:   PetscMPIInt       *N, *displs, l;
1659:   PetscInt           ng, m, i, j, p;
1660:   PetscMPIInt        size;

1662:   PetscFunctionBegin;
1665:   PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)vec), &size));
1666:   PetscCall(VecGetLocalSize(vec, &m));
1667:   PetscCall(VecGetArrayRead(vec, &v));
1668:   PetscCall(PetscMalloc2(m, &tmp, size, &N));
1669:   for (i = 0, j = 0, l = 0; i < m; ++i) {
1670:     /* Can speed this up with sorting */
1671:     for (j = 0; j < l; ++j) {
1672:       if (v[i] == tmp[j]) break;
1673:     }
1674:     if (j == l) {
1675:       tmp[j] = v[i];
1676:       ++l;
1677:     }
1678:   }
1679:   PetscCall(VecRestoreArrayRead(vec, &v));
1680:   /* Gather serial results */
1681:   PetscCallMPI(MPI_Allgather(&l, 1, MPI_INT, N, 1, MPI_INT, PetscObjectComm((PetscObject)vec)));
1682:   for (p = 0, ng = 0; p < size; ++p) ng += N[p];
1683:   PetscCall(PetscMalloc2(ng, &vals, size + 1, &displs));
1684:   for (p = 1, displs[0] = 0; p <= size; ++p) displs[p] = displs[p - 1] + N[p - 1];
1685:   PetscCallMPI(MPI_Allgatherv(tmp, l, MPIU_SCALAR, vals, N, displs, MPIU_SCALAR, PetscObjectComm((PetscObject)vec)));
1686:   /* Find unique entries */
1687: #ifdef PETSC_USE_COMPLEX
1688:   SETERRQ(PetscObjectComm((PetscObject)vec), PETSC_ERR_SUP, "Does not work with complex numbers");
1689: #else
1690:   *n = displs[size];
1691:   PetscCall(PetscSortRemoveDupsReal(n, (PetscReal *)vals));
1692:   if (e) {
1694:     PetscCall(PetscMalloc1(*n, e));
1695:     for (i = 0; i < *n; ++i) (*e)[i] = vals[i];
1696:   }
1697:   PetscCall(PetscFree2(vals, displs));
1698:   PetscCall(PetscFree2(tmp, N));
1699:   PetscFunctionReturn(PETSC_SUCCESS);
1700: #endif
1701: }