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:    Notes:
 19:    One must call VecSetBlockSize() before this routine to set the stride
 20:    information, or use a vector created from a multicomponent DMDA.

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

 24:    Level: advanced

 26: .seealso: `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;

 34:   VecGetLocalSize(v, &n);
 35:   VecGetBlockSize(v, &bs);
 38:   VecGetArray(v, &x);
 39:   for (i = start; i < n; i += bs) x[i] = s;
 40:   VecRestoreArray(v, &x);
 41:   return 0;
 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:    Notes:
 56:    One must call VecSetBlockSize() before this routine to set the stride
 57:    information, or use a vector created from a multicomponent DMDA.

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

 61:    Level: advanced

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

 71:   VecGetLocalSize(v, &n);
 72:   VecGetBlockSize(v, &bs);
 75:   VecGetArray(v, &x);
 76:   for (i = start; i < n; i += bs) x[i] *= scale;
 77:   VecRestoreArray(v, &x);
 78:   return 0;
 79: }

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

 85:    Collective

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

 92:    Output Parameter:
 93: .  norm - the norm

 95:    Notes:
 96:    One must call VecSetBlockSize() before this routine to set the stride
 97:    information, or use a vector created from a multicomponent DMDA.

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

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

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

107:    Level: advanced

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

120:   VecGetLocalSize(v, &n);
121:   VecGetBlockSize(v, &bs);
124:   VecGetArrayRead(v, &x);
125:   if (ntype == NORM_2) {
126:     PetscScalar sum = 0.0;
127:     for (i = start; i < n; i += bs) sum += x[i] * (PetscConj(x[i]));
128:     tnorm = PetscRealPart(sum);
129:     MPIU_Allreduce(&tnorm, nrm, 1, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)v));
130:     *nrm = PetscSqrtReal(*nrm);
131:   } else if (ntype == NORM_1) {
132:     tnorm = 0.0;
133:     for (i = start; i < n; i += bs) tnorm += PetscAbsScalar(x[i]);
134:     MPIU_Allreduce(&tnorm, nrm, 1, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)v));
135:   } else if (ntype == NORM_INFINITY) {
136:     tnorm = 0.0;
137:     for (i = start; i < n; i += bs) {
138:       if (PetscAbsScalar(x[i]) > tnorm) tnorm = PetscAbsScalar(x[i]);
139:     }
140:     MPIU_Allreduce(&tnorm, nrm, 1, MPIU_REAL, MPIU_MAX, PetscObjectComm((PetscObject)v));
141:   } else SETERRQ(PetscObjectComm((PetscObject)v), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown norm type");
142:   VecRestoreArrayRead(v, &x);
143:   return 0;
144: }

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

150:    Collective

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

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

160:    Notes:
161:    One must call VecSetBlockSize() before this routine to set the stride
162:    information, or use a vector created from a multicomponent DMDA.

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

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

171:    Level: advanced

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

183:   VecGetLocalSize(v, &n);
184:   VecGetBlockSize(v, &bs);
187:   VecGetArrayRead(v, &x);
188:   for (i = start; i < n; i += bs) {
189:     if (PetscRealPart(x[i]) > max) {
190:       max = PetscRealPart(x[i]);
191:       id  = i;
192:     }
193:   }
194:   VecRestoreArrayRead(v, &x);
195: #if defined(PETSC_HAVE_MPIUNI)
196:   *nrm = max;
197:   if (idex) *idex = id;
198: #else
199:   if (!idex) {
200:     MPIU_Allreduce(&max, nrm, 1, MPIU_REAL, MPIU_MAX, PetscObjectComm((PetscObject)v));
201:   } else {
202:     struct {
203:       PetscReal v;
204:       PetscInt  i;
205:     } in, out;
206:     PetscInt rstart;

208:     VecGetOwnershipRange(v, &rstart, NULL);
209:     in.v = max;
210:     in.i = rstart + id;
211:     MPIU_Allreduce(&in, &out, 1, MPIU_REAL_INT, MPIU_MAXLOC, PetscObjectComm((PetscObject)v));
212:     *nrm  = out.v;
213:     *idex = out.i;
214:   }
215: #endif
216:   return 0;
217: }

219: /*@
220:    VecStrideMin - Computes the minimum of subvector of a vector defined
221:    by a starting point and a stride and optionally its location.

223:    Collective

225:    Input Parameters:
226: +  v - the vector
227: -  start - starting point of the subvector (defined by a stride)

229:    Output Parameters:
230: +  idex - the location where the minimum occurred. (pass NULL if not required)
231: -  nrm - the minimum value in the subvector

233:    Level: advanced

235:    Notes:
236:    One must call VecSetBlockSize() before this routine to set the stride
237:    information, or use a vector created from a multicomponent DMDA.

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

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

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

256:   VecGetLocalSize(v, &n);
257:   VecGetBlockSize(v, &bs);
260:   VecGetArrayRead(v, &x);
261:   for (i = start; i < n; i += bs) {
262:     if (PetscRealPart(x[i]) < min) {
263:       min = PetscRealPart(x[i]);
264:       id  = i;
265:     }
266:   }
267:   VecRestoreArrayRead(v, &x);
268: #if defined(PETSC_HAVE_MPIUNI)
269:   *nrm = min;
270:   if (idex) *idex = id;
271: #else
272:   if (!idex) {
273:     MPIU_Allreduce(&min, nrm, 1, MPIU_REAL, MPIU_MIN, PetscObjectComm((PetscObject)v));
274:   } else {
275:     struct {
276:       PetscReal v;
277:       PetscInt  i;
278:     } in, out;
279:     PetscInt rstart;

281:     VecGetOwnershipRange(v, &rstart, NULL);
282:     in.v = min;
283:     in.i = rstart + id;
284:     MPIU_Allreduce(&in, &out, 1, MPIU_REAL_INT, MPIU_MINLOC, PetscObjectComm((PetscObject)v));
285:     *nrm  = out.v;
286:     *idex = out.i;
287:   }
288: #endif
289:   return 0;
290: }

292: /*@
293:    VecStrideSum - Computes the sum of subvector of a vector defined
294:    by a starting point and a stride.

296:    Collective

298:    Input Parameters:
299: +  v - the vector
300: -  start - starting point of the subvector (defined by a stride)

302:    Output Parameter:
303: .  sum - the sum

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

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

312:    Level: advanced

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

324:   VecGetLocalSize(v, &n);
325:   VecGetBlockSize(v, &bs);
328:   VecGetArrayRead(v, &x);
329:   for (i = start; i < n; i += bs) local_sum += x[i];
330:   MPIU_Allreduce(&local_sum, sum, 1, MPIU_SCALAR, MPIU_SUM, PetscObjectComm((PetscObject)v));
331:   VecRestoreArrayRead(v, &x);
332:   return 0;
333: }

335: /*@
336:    VecStrideScaleAll - Scales the subvectors of a vector defined
337:    by a starting point and a stride.

339:    Logically Collective

341:    Input Parameters:
342: +  v - the vector
343: -  scales - values to multiply each subvector entry by

345:    Notes:
346:    One must call VecSetBlockSize() before this routine to set the stride
347:    information, or use a vector created from a multicomponent DMDA.

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

351:    Level: advanced

353: .seealso: `VecNorm()`, `VecStrideScale()`, `VecScale()`, `VecStrideGather()`, `VecStrideScatter()`, `VecStrideMin()`, `VecStrideMax()`
354: @*/
355: PetscErrorCode VecStrideScaleAll(Vec v, const PetscScalar *scales)
356: {
357:   PetscInt     i, j, n, bs;
358:   PetscScalar *x;

362:   VecGetLocalSize(v, &n);
363:   VecGetBlockSize(v, &bs);
364:   VecGetArray(v, &x);
365:   /* need to provide optimized code for each bs */
366:   for (i = 0; i < n; i += bs) {
367:     for (j = 0; j < bs; j++) x[i + j] *= scales[j];
368:   }
369:   VecRestoreArray(v, &x);
370:   return 0;
371: }

373: /*@
374:    VecStrideNormAll - Computes the norms of subvectors of a vector defined
375:    by a starting point and a stride.

377:    Collective

379:    Input Parameters:
380: +  v - the vector
381: -  ntype - type of norm, one of NORM_1, NORM_2, NORM_INFINITY

383:    Output Parameter:
384: .  nrm - the norms

386:    Notes:
387:    One must call VecSetBlockSize() before this routine to set the stride
388:    information, or use a vector created from a multicomponent DMDA.

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

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

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

397:    Level: advanced

399: .seealso: `VecNorm()`, `VecStrideGather()`, `VecStrideScatter()`, `VecStrideMin()`, `VecStrideMax()`
400: @*/
401: PetscErrorCode VecStrideNormAll(Vec v, NormType ntype, PetscReal nrm[])
402: {
403:   PetscInt           i, j, n, bs;
404:   const PetscScalar *x;
405:   PetscReal          tnorm[128];
406:   MPI_Comm           comm;

411:   VecGetLocalSize(v, &n);
412:   VecGetArrayRead(v, &x);
413:   PetscObjectGetComm((PetscObject)v, &comm);

415:   VecGetBlockSize(v, &bs);

418:   if (ntype == NORM_2) {
419:     PetscScalar sum[128];
420:     for (j = 0; j < bs; j++) sum[j] = 0.0;
421:     for (i = 0; i < n; i += bs) {
422:       for (j = 0; j < bs; j++) sum[j] += x[i + j] * (PetscConj(x[i + j]));
423:     }
424:     for (j = 0; j < bs; j++) tnorm[j] = PetscRealPart(sum[j]);

426:     MPIU_Allreduce(tnorm, nrm, bs, MPIU_REAL, MPIU_SUM, comm);
427:     for (j = 0; j < bs; j++) nrm[j] = PetscSqrtReal(nrm[j]);
428:   } else if (ntype == NORM_1) {
429:     for (j = 0; j < bs; j++) tnorm[j] = 0.0;

431:     for (i = 0; i < n; i += bs) {
432:       for (j = 0; j < bs; j++) tnorm[j] += PetscAbsScalar(x[i + j]);
433:     }

435:     MPIU_Allreduce(tnorm, nrm, bs, MPIU_REAL, MPIU_SUM, comm);
436:   } else if (ntype == NORM_INFINITY) {
437:     PetscReal tmp;
438:     for (j = 0; j < bs; j++) tnorm[j] = 0.0;

440:     for (i = 0; i < n; i += bs) {
441:       for (j = 0; j < bs; j++) {
442:         if ((tmp = PetscAbsScalar(x[i + j])) > tnorm[j]) tnorm[j] = tmp;
443:         /* check special case of tmp == NaN */
444:         if (tmp != tmp) {
445:           tnorm[j] = tmp;
446:           break;
447:         }
448:       }
449:     }
450:     MPIU_Allreduce(tnorm, nrm, bs, MPIU_REAL, MPIU_MAX, comm);
451:   } else SETERRQ(comm, PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown norm type");
452:   VecRestoreArrayRead(v, &x);
453:   return 0;
454: }

456: /*@
457:    VecStrideMaxAll - Computes the maximums of subvectors of a vector defined
458:    by a starting point and a stride and optionally its location.

460:    Collective

462:    Input Parameter:
463: .  v - the vector

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

470:    Notes:
471:    One must call VecSetBlockSize() before this routine to set the stride
472:    information, or use a vector created from a multicomponent DMDA.

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

476:    Level: advanced

478: .seealso: `VecMax()`, `VecStrideNorm()`, `VecStrideGather()`, `VecStrideScatter()`, `VecStrideMin()`
479: @*/
480: PetscErrorCode VecStrideMaxAll(Vec v, PetscInt idex[], PetscReal nrm[])
481: {
482:   PetscInt           i, j, n, bs;
483:   const PetscScalar *x;
484:   PetscReal          max[128], tmp;
485:   MPI_Comm           comm;

490:   VecGetLocalSize(v, &n);
491:   VecGetArrayRead(v, &x);
492:   PetscObjectGetComm((PetscObject)v, &comm);

494:   VecGetBlockSize(v, &bs);

497:   if (!n) {
498:     for (j = 0; j < bs; j++) max[j] = PETSC_MIN_REAL;
499:   } else {
500:     for (j = 0; j < bs; j++) max[j] = PetscRealPart(x[j]);

502:     for (i = bs; i < n; i += bs) {
503:       for (j = 0; j < bs; j++) {
504:         if ((tmp = PetscRealPart(x[i + j])) > max[j]) max[j] = tmp;
505:       }
506:     }
507:   }
508:   MPIU_Allreduce(max, nrm, bs, MPIU_REAL, MPIU_MAX, comm);

510:   VecRestoreArrayRead(v, &x);
511:   return 0;
512: }

514: /*@
515:    VecStrideMinAll - Computes the minimum of subvector of a vector defined
516:    by a starting point and a stride and optionally its location.

518:    Collective

520:    Input Parameter:
521: .  v - the vector

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

528:    Level: advanced

530:    Notes:
531:    One must call VecSetBlockSize() before this routine to set the stride
532:    information, or use a vector created from a multicomponent DMDA.

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

536: .seealso: `VecMin()`, `VecStrideNorm()`, `VecStrideGather()`, `VecStrideScatter()`, `VecStrideMax()`
537: @*/
538: PetscErrorCode VecStrideMinAll(Vec v, PetscInt idex[], PetscReal nrm[])
539: {
540:   PetscInt           i, n, bs, j;
541:   const PetscScalar *x;
542:   PetscReal          min[128], tmp;
543:   MPI_Comm           comm;

548:   VecGetLocalSize(v, &n);
549:   VecGetArrayRead(v, &x);
550:   PetscObjectGetComm((PetscObject)v, &comm);

552:   VecGetBlockSize(v, &bs);

555:   if (!n) {
556:     for (j = 0; j < bs; j++) min[j] = PETSC_MAX_REAL;
557:   } else {
558:     for (j = 0; j < bs; j++) min[j] = PetscRealPart(x[j]);

560:     for (i = bs; i < n; i += bs) {
561:       for (j = 0; j < bs; j++) {
562:         if ((tmp = PetscRealPart(x[i + j])) < min[j]) min[j] = tmp;
563:       }
564:     }
565:   }
566:   MPIU_Allreduce(min, nrm, bs, MPIU_REAL, MPIU_MIN, comm);

568:   VecRestoreArrayRead(v, &x);
569:   return 0;
570: }

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

575:    Collective

577:    Input Parameter:
578: .  v - the vector

580:    Output Parameter:
581: .  sums - the sums

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

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

590:    Level: advanced

592: .seealso: `VecSum()`, `VecStrideGather()`, `VecStrideScatter()`, `VecStrideMin()`, `VecStrideMax()`
593: @*/
594: PetscErrorCode VecStrideSumAll(Vec v, PetscScalar sums[])
595: {
596:   PetscInt           i, j, n, bs;
597:   const PetscScalar *x;
598:   PetscScalar        local_sums[128];
599:   MPI_Comm           comm;

603:   VecGetLocalSize(v, &n);
604:   VecGetArrayRead(v, &x);
605:   PetscObjectGetComm((PetscObject)v, &comm);

607:   VecGetBlockSize(v, &bs);

610:   for (j = 0; j < bs; j++) local_sums[j] = 0.0;
611:   for (i = 0; i < n; i += bs) {
612:     for (j = 0; j < bs; j++) local_sums[j] += x[i + j];
613:   }
614:   MPIU_Allreduce(local_sums, sums, bs, MPIU_SCALAR, MPIU_SUM, comm);

616:   VecRestoreArrayRead(v, &x);
617:   return 0;
618: }

620: /*----------------------------------------------------------------------------------------------*/
621: /*@
622:    VecStrideGatherAll - Gathers all the single components from a multi-component vector into
623:    separate vectors.

625:    Collective

627:    Input Parameters:
628: +  v - the vector
629: -  addv - one of ADD_VALUES,INSERT_VALUES,MAX_VALUES

631:    Output Parameter:
632: .  s - the location where the subvectors are stored

634:    Notes:
635:    One must call VecSetBlockSize() before this routine to set the stride
636:    information, or use a vector created from a multicomponent DMDA.

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

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

645:    Not optimized; could be easily

647:    Level: advanced

649: .seealso: `VecStrideNorm()`, `VecStrideScatter()`, `VecStrideMin()`, `VecStrideMax()`, `VecStrideGather()`,
650:           `VecStrideScatterAll()`
651: @*/
652: PetscErrorCode VecStrideGatherAll(Vec v, Vec s[], InsertMode addv)
653: {
654:   PetscInt           i, n, n2, bs, j, k, *bss = NULL, nv, jj, nvc;
655:   PetscScalar      **y;
656:   const PetscScalar *x;

661:   VecGetLocalSize(v, &n);
662:   VecGetLocalSize(s[0], &n2);
663:   VecGetArrayRead(v, &x);
664:   VecGetBlockSize(v, &bs);

667:   PetscMalloc2(bs, &y, bs, &bss);
668:   nv  = 0;
669:   nvc = 0;
670:   for (i = 0; i < bs; i++) {
671:     VecGetBlockSize(s[i], &bss[i]);
672:     if (bss[i] < 1) bss[i] = 1; /* if user never set it then assume 1  Re: [PETSC #8241] VecStrideGatherAll */
673:     VecGetArray(s[i], &y[i]);
674:     nvc += bss[i];
675:     nv++;
677:     if (nvc == bs) break;
678:   }

680:   n = n / bs;

682:   jj = 0;
683:   if (addv == INSERT_VALUES) {
684:     for (j = 0; j < nv; j++) {
685:       for (k = 0; k < bss[j]; k++) {
686:         for (i = 0; i < n; i++) y[j][i * bss[j] + k] = x[bs * i + jj + k];
687:       }
688:       jj += bss[j];
689:     }
690:   } else if (addv == ADD_VALUES) {
691:     for (j = 0; j < nv; j++) {
692:       for (k = 0; k < bss[j]; k++) {
693:         for (i = 0; i < n; i++) y[j][i * bss[j] + k] += x[bs * i + jj + k];
694:       }
695:       jj += bss[j];
696:     }
697: #if !defined(PETSC_USE_COMPLEX)
698:   } else if (addv == MAX_VALUES) {
699:     for (j = 0; j < nv; j++) {
700:       for (k = 0; k < bss[j]; k++) {
701:         for (i = 0; i < n; i++) y[j][i * bss[j] + k] = PetscMax(y[j][i * bss[j] + k], x[bs * i + jj + k]);
702:       }
703:       jj += bss[j];
704:     }
705: #endif
706:   } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown insert type");

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

711:   PetscFree2(y, bss);
712:   return 0;
713: }

715: /*@
716:    VecStrideScatterAll - Scatters all the single components from separate vectors into
717:      a multi-component vector.

719:    Collective

721:    Input Parameters:
722: +  s - the location where the subvectors are stored
723: -  addv - one of ADD_VALUES,INSERT_VALUES,MAX_VALUES

725:    Output Parameter:
726: .  v - the multicomponent vector

728:    Notes:
729:    One must call VecSetBlockSize() before this routine to set the stride
730:    information, or use a vector created from a multicomponent DMDA.

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

735:    Not optimized; could be easily

737:    Level: advanced

739: .seealso: `VecStrideNorm()`, `VecStrideScatter()`, `VecStrideMin()`, `VecStrideMax()`, `VecStrideGather()`,
740:           `VecStrideScatterAll()`
741: @*/
742: PetscErrorCode VecStrideScatterAll(Vec s[], Vec v, InsertMode addv)
743: {
744:   PetscInt            i, n, n2, bs, j, jj, k, *bss = NULL, nv, nvc;
745:   PetscScalar        *x;
746:   PetscScalar const **y;

751:   VecGetLocalSize(v, &n);
752:   VecGetLocalSize(s[0], &n2);
753:   VecGetArray(v, &x);
754:   VecGetBlockSize(v, &bs);

757:   PetscMalloc2(bs, (PetscScalar ***)&y, bs, &bss);
758:   nv  = 0;
759:   nvc = 0;
760:   for (i = 0; i < bs; i++) {
761:     VecGetBlockSize(s[i], &bss[i]);
762:     if (bss[i] < 1) bss[i] = 1; /* if user never set it then assume 1  Re: [PETSC #8241] VecStrideGatherAll */
763:     VecGetArrayRead(s[i], &y[i]);
764:     nvc += bss[i];
765:     nv++;
767:     if (nvc == bs) break;
768:   }

770:   n = n / bs;

772:   jj = 0;
773:   if (addv == INSERT_VALUES) {
774:     for (j = 0; j < nv; j++) {
775:       for (k = 0; k < bss[j]; k++) {
776:         for (i = 0; i < n; i++) x[bs * i + jj + k] = y[j][i * bss[j] + k];
777:       }
778:       jj += bss[j];
779:     }
780:   } else if (addv == ADD_VALUES) {
781:     for (j = 0; j < nv; j++) {
782:       for (k = 0; k < bss[j]; k++) {
783:         for (i = 0; i < n; i++) x[bs * i + jj + k] += y[j][i * bss[j] + k];
784:       }
785:       jj += bss[j];
786:     }
787: #if !defined(PETSC_USE_COMPLEX)
788:   } else if (addv == MAX_VALUES) {
789:     for (j = 0; j < nv; j++) {
790:       for (k = 0; k < bss[j]; k++) {
791:         for (i = 0; i < n; i++) x[bs * i + jj + k] = PetscMax(x[bs * i + jj + k], y[j][i * bss[j] + k]);
792:       }
793:       jj += bss[j];
794:     }
795: #endif
796:   } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown insert type");

798:   VecRestoreArray(v, &x);
799:   for (i = 0; i < nv; i++) VecRestoreArrayRead(s[i], &y[i]);
800:   PetscFree2(*(PetscScalar ***)&y, bss);
801:   return 0;
802: }

804: /*@
805:    VecStrideGather - Gathers a single component from a multi-component vector into
806:    another vector.

808:    Collective

810:    Input Parameters:
811: +  v - the vector
812: .  start - starting point of the subvector (defined by a stride)
813: -  addv - one of ADD_VALUES,INSERT_VALUES,MAX_VALUES

815:    Output Parameter:
816: .  s - the location where the subvector is stored

818:    Notes:
819:    One must call VecSetBlockSize() before this routine to set the stride
820:    information, or use a vector created from a multicomponent DMDA.

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

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

828:    Not optimized; could be easily

830:    Level: advanced

832: .seealso: `VecStrideNorm()`, `VecStrideScatter()`, `VecStrideMin()`, `VecStrideMax()`, `VecStrideGatherAll()`,
833:           `VecStrideScatterAll()`
834: @*/
835: PetscErrorCode VecStrideGather(Vec v, PetscInt start, Vec s, InsertMode addv)
836: {
841:              v->map->bs);
842:   PetscUseTypeMethod(v, stridegather, start, s, addv);
843:   return 0;
844: }

846: /*@
847:    VecStrideScatter - Scatters a single component from a vector into a multi-component vector.

849:    Collective

851:    Input Parameters:
852: +  s - the single-component vector
853: .  start - starting point of the subvector (defined by a stride)
854: -  addv - one of ADD_VALUES,INSERT_VALUES,MAX_VALUES

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

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

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

866:    Not optimized; could be easily

868:    Level: advanced

870: .seealso: `VecStrideNorm()`, `VecStrideGather()`, `VecStrideMin()`, `VecStrideMax()`, `VecStrideGatherAll()`,
871:           `VecStrideScatterAll()`, `VecStrideSubSetScatter()`, `VecStrideSubSetGather()`
872: @*/
873: PetscErrorCode VecStrideScatter(Vec s, PetscInt start, Vec v, InsertMode addv)
874: {
879:              v->map->bs);
880:   (*v->ops->stridescatter)(s, start, v, addv);
881:   return 0;
882: }

884: /*@
885:    VecStrideSubSetGather - Gathers a subset of components from a multi-component vector into
886:    another vector.

888:    Collective

890:    Input Parameters:
891: +  v - the vector
892: .  nidx - the number of indices
893: .  idxv - the indices of the components 0 <= idxv[0] ...idxv[nidx-1] < bs(v), they need not be sorted
894: .  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
895: -  addv - one of ADD_VALUES,INSERT_VALUES,MAX_VALUES

897:    Output Parameter:
898: .  s - the location where the subvector is stored

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

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

906:    Not optimized; could be easily

908:    Level: advanced

910: .seealso: `VecStrideNorm()`, `VecStrideScatter()`, `VecStrideGather()`, `VecStrideSubSetScatter()`, `VecStrideMin()`, `VecStrideMax()`, `VecStrideGatherAll()`,
911:           `VecStrideScatterAll()`
912: @*/
913: PetscErrorCode VecStrideSubSetGather(Vec v, PetscInt nidx, const PetscInt idxv[], const PetscInt idxs[], Vec s, InsertMode addv)
914: {
917:   if (nidx == PETSC_DETERMINE) nidx = s->map->bs;
918:   PetscUseTypeMethod(v, stridesubsetgather, nidx, idxv, idxs, s, addv);
919:   return 0;
920: }

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

925:    Collective

927:    Input Parameters:
928: +  s - the smaller-component vector
929: .  nidx - the number of indices in idx
930: .  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
931: .  idxv - the indices of the components in the larger-component vector, 0 <= idx[0] ...idx[nidx-1] < bs(v) they need not be sorted
932: -  addv - one of ADD_VALUES,INSERT_VALUES,MAX_VALUES

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

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

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

943:    Not optimized; could be easily

945:    Level: advanced

947: .seealso: `VecStrideNorm()`, `VecStrideGather()`, `VecStrideGather()`, `VecStrideSubSetGather()`, `VecStrideMin()`, `VecStrideMax()`, `VecStrideGatherAll()`,
948:           `VecStrideScatterAll()`
949: @*/
950: PetscErrorCode VecStrideSubSetScatter(Vec s, PetscInt nidx, const PetscInt idxs[], const PetscInt idxv[], Vec v, InsertMode addv)
951: {
954:   if (nidx == PETSC_DETERMINE) nidx = s->map->bs;
955:   (*v->ops->stridesubsetscatter)(s, nidx, idxs, idxv, v, addv);
956:   return 0;
957: }

959: PetscErrorCode VecStrideGather_Default(Vec v, PetscInt start, Vec s, InsertMode addv)
960: {
961:   PetscInt           i, n, bs, ns;
962:   const PetscScalar *x;
963:   PetscScalar       *y;

965:   VecGetLocalSize(v, &n);
966:   VecGetLocalSize(s, &ns);
967:   VecGetArrayRead(v, &x);
968:   VecGetArray(s, &y);

970:   bs = v->map->bs;
972:   x += start;
973:   n = n / bs;

975:   if (addv == INSERT_VALUES) {
976:     for (i = 0; i < n; i++) y[i] = x[bs * i];
977:   } else if (addv == ADD_VALUES) {
978:     for (i = 0; i < n; i++) y[i] += x[bs * i];
979: #if !defined(PETSC_USE_COMPLEX)
980:   } else if (addv == MAX_VALUES) {
981:     for (i = 0; i < n; i++) y[i] = PetscMax(y[i], x[bs * i]);
982: #endif
983:   } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown insert type");

985:   VecRestoreArrayRead(v, &x);
986:   VecRestoreArray(s, &y);
987:   return 0;
988: }

990: PetscErrorCode VecStrideScatter_Default(Vec s, PetscInt start, Vec v, InsertMode addv)
991: {
992:   PetscInt           i, n, bs, ns;
993:   PetscScalar       *x;
994:   const PetscScalar *y;

996:   VecGetLocalSize(v, &n);
997:   VecGetLocalSize(s, &ns);
998:   VecGetArray(v, &x);
999:   VecGetArrayRead(s, &y);

1001:   bs = v->map->bs;
1003:   x += start;
1004:   n = n / bs;

1006:   if (addv == INSERT_VALUES) {
1007:     for (i = 0; i < n; i++) x[bs * i] = y[i];
1008:   } else if (addv == ADD_VALUES) {
1009:     for (i = 0; i < n; i++) x[bs * i] += y[i];
1010: #if !defined(PETSC_USE_COMPLEX)
1011:   } else if (addv == MAX_VALUES) {
1012:     for (i = 0; i < n; i++) x[bs * i] = PetscMax(y[i], x[bs * i]);
1013: #endif
1014:   } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown insert type");

1016:   VecRestoreArray(v, &x);
1017:   VecRestoreArrayRead(s, &y);
1018:   return 0;
1019: }

1021: PetscErrorCode VecStrideSubSetGather_Default(Vec v, PetscInt nidx, const PetscInt idxv[], const PetscInt idxs[], Vec s, InsertMode addv)
1022: {
1023:   PetscInt           i, j, n, bs, bss, ns;
1024:   const PetscScalar *x;
1025:   PetscScalar       *y;

1027:   VecGetLocalSize(v, &n);
1028:   VecGetLocalSize(s, &ns);
1029:   VecGetArrayRead(v, &x);
1030:   VecGetArray(s, &y);

1032:   bs  = v->map->bs;
1033:   bss = s->map->bs;
1034:   n   = n / bs;

1036:   if (PetscDefined(USE_DEBUG)) {
1038:     for (j = 0; j < nidx; j++) {
1041:     }
1043:   }

1045:   if (addv == INSERT_VALUES) {
1046:     if (!idxs) {
1047:       for (i = 0; i < n; i++) {
1048:         for (j = 0; j < bss; j++) y[bss * i + j] = x[bs * i + idxv[j]];
1049:       }
1050:     } else {
1051:       for (i = 0; i < n; i++) {
1052:         for (j = 0; j < bss; j++) y[bss * i + idxs[j]] = x[bs * i + idxv[j]];
1053:       }
1054:     }
1055:   } else if (addv == ADD_VALUES) {
1056:     if (!idxs) {
1057:       for (i = 0; i < n; i++) {
1058:         for (j = 0; j < bss; j++) y[bss * i + j] += x[bs * i + idxv[j]];
1059:       }
1060:     } else {
1061:       for (i = 0; i < n; i++) {
1062:         for (j = 0; j < bss; j++) y[bss * i + idxs[j]] += x[bs * i + idxv[j]];
1063:       }
1064:     }
1065: #if !defined(PETSC_USE_COMPLEX)
1066:   } else if (addv == MAX_VALUES) {
1067:     if (!idxs) {
1068:       for (i = 0; i < n; i++) {
1069:         for (j = 0; j < bss; j++) y[bss * i + j] = PetscMax(y[bss * i + j], x[bs * i + idxv[j]]);
1070:       }
1071:     } else {
1072:       for (i = 0; i < n; i++) {
1073:         for (j = 0; j < bss; j++) y[bss * i + idxs[j]] = PetscMax(y[bss * i + idxs[j]], x[bs * i + idxv[j]]);
1074:       }
1075:     }
1076: #endif
1077:   } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown insert type");

1079:   VecRestoreArrayRead(v, &x);
1080:   VecRestoreArray(s, &y);
1081:   return 0;
1082: }

1084: PetscErrorCode VecStrideSubSetScatter_Default(Vec s, PetscInt nidx, const PetscInt idxs[], const PetscInt idxv[], Vec v, InsertMode addv)
1085: {
1086:   PetscInt           j, i, n, bs, ns, bss;
1087:   PetscScalar       *x;
1088:   const PetscScalar *y;

1090:   VecGetLocalSize(v, &n);
1091:   VecGetLocalSize(s, &ns);
1092:   VecGetArray(v, &x);
1093:   VecGetArrayRead(s, &y);

1095:   bs  = v->map->bs;
1096:   bss = s->map->bs;
1097:   n   = n / bs;

1099:   if (PetscDefined(USE_DEBUG)) {
1101:     for (j = 0; j < bss; j++) {
1102:       if (idxs) {
1105:       }
1106:     }
1108:   }

1110:   if (addv == INSERT_VALUES) {
1111:     if (!idxs) {
1112:       for (i = 0; i < n; i++) {
1113:         for (j = 0; j < bss; j++) x[bs * i + idxv[j]] = y[bss * i + j];
1114:       }
1115:     } else {
1116:       for (i = 0; i < n; i++) {
1117:         for (j = 0; j < bss; j++) x[bs * i + idxv[j]] = y[bss * i + idxs[j]];
1118:       }
1119:     }
1120:   } else if (addv == ADD_VALUES) {
1121:     if (!idxs) {
1122:       for (i = 0; i < n; i++) {
1123:         for (j = 0; j < bss; j++) x[bs * i + idxv[j]] += y[bss * i + j];
1124:       }
1125:     } else {
1126:       for (i = 0; i < n; i++) {
1127:         for (j = 0; j < bss; j++) x[bs * i + idxv[j]] += y[bss * i + idxs[j]];
1128:       }
1129:     }
1130: #if !defined(PETSC_USE_COMPLEX)
1131:   } else if (addv == MAX_VALUES) {
1132:     if (!idxs) {
1133:       for (i = 0; i < n; i++) {
1134:         for (j = 0; j < bss; j++) x[bs * i + idxv[j]] = PetscMax(y[bss * i + j], x[bs * i + idxv[j]]);
1135:       }
1136:     } else {
1137:       for (i = 0; i < n; i++) {
1138:         for (j = 0; j < bss; j++) x[bs * i + idxv[j]] = PetscMax(y[bss * i + idxs[j]], x[bs * i + idxv[j]]);
1139:       }
1140:     }
1141: #endif
1142:   } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown insert type");

1144:   VecRestoreArray(v, &x);
1145:   VecRestoreArrayRead(s, &y);
1146:   return 0;
1147: }

1149: PetscErrorCode VecReciprocal_Default(Vec v)
1150: {
1151:   PetscInt     i, n;
1152:   PetscScalar *x;

1154:   VecGetLocalSize(v, &n);
1155:   VecGetArray(v, &x);
1156:   for (i = 0; i < n; i++) {
1157:     if (x[i] != (PetscScalar)0.0) x[i] = (PetscScalar)1.0 / x[i];
1158:   }
1159:   VecRestoreArray(v, &x);
1160:   return 0;
1161: }

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

1166:   Not collective

1168:   Input Parameter:
1169: . v - The vector

1171:   Output Parameter:
1172: . v - The vector of exponents

1174:   Level: beginner

1176: .seealso: `VecLog()`, `VecAbs()`, `VecSqrtAbs()`, `VecReciprocal()`

1178: @*/
1179: PetscErrorCode VecExp(Vec v)
1180: {
1181:   PetscScalar *x;
1182:   PetscInt     i, n;

1185:   if (v->ops->exp) PetscUseTypeMethod(v, exp);
1186:   else {
1187:     VecGetLocalSize(v, &n);
1188:     VecGetArray(v, &x);
1189:     for (i = 0; i < n; i++) x[i] = PetscExpScalar(x[i]);
1190:     VecRestoreArray(v, &x);
1191:   }
1192:   return 0;
1193: }

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

1198:   Not collective

1200:   Input Parameter:
1201: . v - The vector

1203:   Output Parameter:
1204: . v - The vector of logs

1206:   Level: beginner

1208: .seealso: `VecExp()`, `VecAbs()`, `VecSqrtAbs()`, `VecReciprocal()`

1210: @*/
1211: PetscErrorCode VecLog(Vec v)
1212: {
1213:   PetscScalar *x;
1214:   PetscInt     i, n;

1217:   if (v->ops->log) PetscUseTypeMethod(v, log);
1218:   else {
1219:     VecGetLocalSize(v, &n);
1220:     VecGetArray(v, &x);
1221:     for (i = 0; i < n; i++) x[i] = PetscLogScalar(x[i]);
1222:     VecRestoreArray(v, &x);
1223:   }
1224:   return 0;
1225: }

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

1230:   Not collective

1232:   Input Parameter:
1233: . v - The vector

1235:   Output Parameter:
1236: . v - The vector square root

1238:   Level: beginner

1240:   Note: The actual function is sqrt(|x_i|)

1242: .seealso: `VecLog()`, `VecExp()`, `VecReciprocal()`, `VecAbs()`

1244: @*/
1245: PetscErrorCode VecSqrtAbs(Vec v)
1246: {
1247:   PetscScalar *x;
1248:   PetscInt     i, n;

1251:   if (v->ops->sqrt) PetscUseTypeMethod(v, sqrt);
1252:   else {
1253:     VecGetLocalSize(v, &n);
1254:     VecGetArray(v, &x);
1255:     for (i = 0; i < n; i++) x[i] = PetscSqrtReal(PetscAbsScalar(x[i]));
1256:     VecRestoreArray(v, &x);
1257:   }
1258:   return 0;
1259: }

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

1264:   Collective

1266:   Input Parameters:
1267: + s - first vector
1268: - t - second vector

1270:   Output Parameters:
1271: + dp - s'conj(t)
1272: - nm - t'conj(t)

1274:   Level: advanced

1276:   Notes:
1277:     conj(x) is the complex conjugate of x when x is complex

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

1281: @*/
1282: PetscErrorCode VecDotNorm2(Vec s, Vec t, PetscScalar *dp, PetscReal *nm)
1283: {
1284:   const PetscScalar *sx, *tx;
1285:   PetscScalar        dpx = 0.0, nmx = 0.0, work[2], sum[2];
1286:   PetscInt           i, n;


1298:   PetscLogEventBegin(VEC_DotNorm2, s, t, 0, 0);
1299:   if (s->ops->dotnorm2) {
1300:     PetscUseTypeMethod(s, dotnorm2, t, dp, &dpx);
1301:     *nm = PetscRealPart(dpx);
1302:   } else {
1303:     VecGetLocalSize(s, &n);
1304:     VecGetArrayRead(s, &sx);
1305:     VecGetArrayRead(t, &tx);

1307:     for (i = 0; i < n; i++) {
1308:       dpx += sx[i] * PetscConj(tx[i]);
1309:       nmx += tx[i] * PetscConj(tx[i]);
1310:     }
1311:     work[0] = dpx;
1312:     work[1] = nmx;

1314:     MPIU_Allreduce(work, sum, 2, MPIU_SCALAR, MPIU_SUM, PetscObjectComm((PetscObject)s));
1315:     *dp = sum[0];
1316:     *nm = PetscRealPart(sum[1]);

1318:     VecRestoreArrayRead(t, &tx);
1319:     VecRestoreArrayRead(s, &sx);
1320:     PetscLogFlops(4.0 * n);
1321:   }
1322:   PetscLogEventEnd(VEC_DotNorm2, s, t, 0, 0);
1323:   return 0;
1324: }

1326: /*@
1327:    VecSum - Computes the sum of all the components of a vector.

1329:    Collective

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

1334:    Output Parameter:
1335: .  sum - the result

1337:    Level: beginner

1339: .seealso: `VecMean()`, `VecNorm()`
1340: @*/
1341: PetscErrorCode VecSum(Vec v, PetscScalar *sum)
1342: {
1343:   PetscInt           i, n;
1344:   const PetscScalar *x;

1348:   *sum = 0.0;
1349:   if (v->ops->sum) {
1350:     PetscUseTypeMethod(v, sum, sum);
1351:   } else {
1352:     VecGetLocalSize(v, &n);
1353:     VecGetArrayRead(v, &x);
1354:     for (i = 0; i < n; i++) *sum += x[i];
1355:     VecRestoreArrayRead(v, &x);
1356:   }
1357:   MPIU_Allreduce(MPI_IN_PLACE, sum, 1, MPIU_SCALAR, MPIU_SUM, PetscObjectComm((PetscObject)v));
1358:   return 0;
1359: }

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

1364:    Collective

1366:    Input Parameter:
1367: .  v - the vector

1369:    Output Parameter:
1370: .  mean - the result

1372:    Level: beginner

1374: .seealso: `VecSum()`, `VecNorm()`
1375: @*/
1376: PetscErrorCode VecMean(Vec v, PetscScalar *mean)
1377: {
1378:   PetscInt n;

1382:   VecGetSize(v, &n);
1383:   VecSum(v, mean);
1384:   *mean /= n;
1385:   return 0;
1386: }

1388: /*@
1389:    VecImaginaryPart - Replaces a complex vector with its imginary part

1391:    Collective

1393:    Input Parameter:
1394: .  v - the vector

1396:    Level: beginner

1398: .seealso: `VecNorm()`, `VecRealPart()`
1399: @*/
1400: PetscErrorCode VecImaginaryPart(Vec v)
1401: {
1402:   PetscInt     i, n;
1403:   PetscScalar *x;

1406:   VecGetLocalSize(v, &n);
1407:   VecGetArray(v, &x);
1408:   for (i = 0; i < n; i++) x[i] = PetscImaginaryPart(x[i]);
1409:   VecRestoreArray(v, &x);
1410:   return 0;
1411: }

1413: /*@
1414:    VecRealPart - Replaces a complex vector with its real part

1416:    Collective

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

1421:    Level: beginner

1423: .seealso: `VecNorm()`, `VecImaginaryPart()`
1424: @*/
1425: PetscErrorCode VecRealPart(Vec v)
1426: {
1427:   PetscInt     i, n;
1428:   PetscScalar *x;

1431:   VecGetLocalSize(v, &n);
1432:   VecGetArray(v, &x);
1433:   for (i = 0; i < n; i++) x[i] = PetscRealPart(x[i]);
1434:   VecRestoreArray(v, &x);
1435:   return 0;
1436: }

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

1442:    Logically Collective

1444:    Input Parameters:
1445: +  v - the vector
1446: -  shift - the shift

1448:    Level: intermediate

1450: @*/
1451: PetscErrorCode VecShift(Vec v, PetscScalar shift)
1452: {
1453:   PetscInt     i, n;
1454:   PetscScalar *x;

1458:   VecSetErrorIfLocked(v, 1);
1459:   if (shift == 0.0) return 0;

1461:   if (v->ops->shift) PetscUseTypeMethod(v, shift, shift);
1462:   else {
1463:     VecGetLocalSize(v, &n);
1464:     VecGetArray(v, &x);
1465:     for (i = 0; i < n; i++) x[i] += shift;
1466:     VecRestoreArray(v, &x);
1467:   }
1468:   return 0;
1469: }

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

1474:    Logically Collective

1476:    Input Parameters:
1477: .  v - the vector

1479:    Level: intermediate

1481: @*/
1482: PetscErrorCode VecAbs(Vec v)
1483: {
1484:   PetscInt     i, n;
1485:   PetscScalar *x;

1488:   VecSetErrorIfLocked(v, 1);

1490:   if (v->ops->abs) {
1491:     PetscUseTypeMethod(v, abs);
1492:   } else {
1493:     VecGetLocalSize(v, &n);
1494:     VecGetArray(v, &x);
1495:     for (i = 0; i < n; i++) x[i] = PetscAbsScalar(x[i]);
1496:     VecRestoreArray(v, &x);
1497:   }
1498:   return 0;
1499: }

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

1504:   Input Parameters:
1505: + vec   - The vector
1506: . order - The ordering
1507: - inv   - The flag for inverting the permutation

1509:   Level: beginner

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

1513: .seealso: `MatPermute()`
1514: @*/
1515: PetscErrorCode VecPermute(Vec x, IS row, PetscBool inv)
1516: {
1517:   const PetscScalar *array;
1518:   PetscScalar       *newArray;
1519:   const PetscInt    *idx;
1520:   PetscInt           i, rstart, rend;

1524:   VecSetErrorIfLocked(x, 1);
1525:   VecGetOwnershipRange(x, &rstart, &rend);
1526:   ISGetIndices(row, &idx);
1527:   VecGetArrayRead(x, &array);
1528:   PetscMalloc1(x->map->n, &newArray);
1529:   if (PetscDefined(USE_DEBUG)) {
1531:   }
1532:   if (!inv) {
1533:     for (i = 0; i < x->map->n; i++) newArray[i] = array[idx[i] - rstart];
1534:   } else {
1535:     for (i = 0; i < x->map->n; i++) newArray[idx[i] - rstart] = array[i];
1536:   }
1537:   VecRestoreArrayRead(x, &array);
1538:   ISRestoreIndices(row, &idx);
1539:   VecReplaceArray(x, newArray);
1540:   return 0;
1541: }

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

1548:    Collective

1550:    Input Parameters:
1551: +  vec1 - the first vector
1552: -  vec2 - the second vector

1554:    Output Parameter:
1555: .  flg - PETSC_TRUE if the vectors are equal; PETSC_FALSE otherwise.

1557:    Level: intermediate
1558: @*/
1559: PetscErrorCode VecEqual(Vec vec1, Vec vec2, PetscBool *flg)
1560: {
1561:   const PetscScalar *v1, *v2;
1562:   PetscInt           n1, n2, N1, N2;
1563:   PetscBool          flg1;

1568:   if (vec1 == vec2) *flg = PETSC_TRUE;
1569:   else {
1570:     VecGetSize(vec1, &N1);
1571:     VecGetSize(vec2, &N2);
1572:     if (N1 != N2) flg1 = PETSC_FALSE;
1573:     else {
1574:       VecGetLocalSize(vec1, &n1);
1575:       VecGetLocalSize(vec2, &n2);
1576:       if (n1 != n2) flg1 = PETSC_FALSE;
1577:       else {
1578:         VecGetArrayRead(vec1, &v1);
1579:         VecGetArrayRead(vec2, &v2);
1580:         PetscArraycmp(v1, v2, n1, &flg1);
1581:         VecRestoreArrayRead(vec1, &v1);
1582:         VecRestoreArrayRead(vec2, &v2);
1583:       }
1584:     }
1585:     /* combine results from all processors */
1586:     MPIU_Allreduce(&flg1, flg, 1, MPIU_BOOL, MPI_MIN, PetscObjectComm((PetscObject)vec1));
1587:   }
1588:   return 0;
1589: }

1591: /*@
1592:    VecUniqueEntries - Compute the number of unique entries, and those entries

1594:    Collective

1596:    Input Parameter:
1597: .  vec - the vector

1599:    Output Parameters:
1600: +  n - The number of unique entries
1601: -  e - The entries

1603:    Level: intermediate

1605: @*/
1606: PetscErrorCode VecUniqueEntries(Vec vec, PetscInt *n, PetscScalar **e)
1607: {
1608:   const PetscScalar *v;
1609:   PetscScalar       *tmp, *vals;
1610:   PetscMPIInt       *N, *displs, l;
1611:   PetscInt           ng, m, i, j, p;
1612:   PetscMPIInt        size;

1616:   MPI_Comm_size(PetscObjectComm((PetscObject)vec), &size);
1617:   VecGetLocalSize(vec, &m);
1618:   VecGetArrayRead(vec, &v);
1619:   PetscMalloc2(m, &tmp, size, &N);
1620:   for (i = 0, j = 0, l = 0; i < m; ++i) {
1621:     /* Can speed this up with sorting */
1622:     for (j = 0; j < l; ++j) {
1623:       if (v[i] == tmp[j]) break;
1624:     }
1625:     if (j == l) {
1626:       tmp[j] = v[i];
1627:       ++l;
1628:     }
1629:   }
1630:   VecRestoreArrayRead(vec, &v);
1631:   /* Gather serial results */
1632:   MPI_Allgather(&l, 1, MPI_INT, N, 1, MPI_INT, PetscObjectComm((PetscObject)vec));
1633:   for (p = 0, ng = 0; p < size; ++p) ng += N[p];
1634:   PetscMalloc2(ng, &vals, size + 1, &displs);
1635:   for (p = 1, displs[0] = 0; p <= size; ++p) displs[p] = displs[p - 1] + N[p - 1];
1636:   MPI_Allgatherv(tmp, l, MPIU_SCALAR, vals, N, displs, MPIU_SCALAR, PetscObjectComm((PetscObject)vec));
1637:   /* Find unique entries */
1638: #ifdef PETSC_USE_COMPLEX
1639:   SETERRQ(PetscObjectComm((PetscObject)vec), PETSC_ERR_SUP, "Does not work with complex numbers");
1640: #else
1641:   *n = displs[size];
1642:   PetscSortRemoveDupsReal(n, (PetscReal *)vals);
1643:   if (e) {
1645:     PetscMalloc1(*n, e);
1646:     for (i = 0; i < *n; ++i) (*e)[i] = vals[i];
1647:   }
1648:   PetscFree2(vals, displs);
1649:   PetscFree2(tmp, N);
1650:   return 0;
1651: #endif
1652: }