Actual source code: vinv.c

petsc-main 2021-04-15
Report Typos and Errors

  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 on Vec

 13:    Input Parameter:
 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


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

 39:   VecSetErrorIfLocked(v,1);

 41:   VecGetLocalSize(v,&n);
 42:   VecGetArray(v,&x);
 43:   VecGetBlockSize(v,&bs);
 44:   if (start < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative start %D",start);
 45:   else if (start >= bs) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Start of stride subvector (%D) is too large for stride\n  Have you set the vector blocksize (%D) correctly with VecSetBlockSize()?",start,bs);
 46:   x += start;
 47:   for (i=0; i<n; i+=bs) x[i] = s;
 48:   x -= start;
 49:   VecRestoreArray(v,&x);
 50:   return(0);
 51: }

 53: /*@
 54:    VecStrideScale - Scales a subvector of a vector defined
 55:    by a starting point and a stride.

 57:    Logically Collective on Vec

 59:    Input Parameter:
 60: +  v - the vector
 61: .  start - starting point of the subvector (defined by a stride)
 62: -  scale - value to multiply each subvector entry by

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

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

 70:    Level: advanced


 73: .seealso: VecNorm(), VecStrideGather(), VecStrideScatter(), VecStrideMin(), VecStrideMax(), VecStrideScale()
 74: @*/
 75: PetscErrorCode  VecStrideScale(Vec v,PetscInt start,PetscScalar scale)
 76: {
 78:   PetscInt       i,n,bs;
 79:   PetscScalar    *x;

 85:   VecSetErrorIfLocked(v,1);

 87:   VecGetLocalSize(v,&n);
 88:   VecGetArray(v,&x);
 89:   VecGetBlockSize(v,&bs);
 90:   if (start < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative start %D",start);
 91:   else if (start >= bs) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Start of stride subvector (%D) is too large for stride\n  Have you set the vector blocksize (%D) correctly with VecSetBlockSize()?",start,bs);
 92:   x += start;
 93:   for (i=0; i<n; i+=bs) x[i] *= scale;
 94:   x -= start;
 95:   VecRestoreArray(v,&x);
 96:   return(0);
 97: }

 99: /*@
100:    VecStrideNorm - Computes the norm of subvector of a vector defined
101:    by a starting point and a stride.

103:    Collective on Vec

105:    Input Parameter:
106: +  v - the vector
107: .  start - starting point of the subvector (defined by a stride)
108: -  ntype - type of norm, one of NORM_1, NORM_2, NORM_INFINITY

110:    Output Parameter:
111: .  norm - the norm

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

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

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

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

125:    Level: advanced


128: .seealso: VecNorm(), VecStrideGather(), VecStrideScatter(), VecStrideMin(), VecStrideMax()
129: @*/
130: PetscErrorCode  VecStrideNorm(Vec v,PetscInt start,NormType ntype,PetscReal *nrm)
131: {
132:   PetscErrorCode    ierr;
133:   PetscInt          i,n,bs;
134:   const PetscScalar *x;
135:   PetscReal         tnorm;
136:   MPI_Comm          comm;

141:   VecGetLocalSize(v,&n);
142:   VecGetArrayRead(v,&x);
143:   PetscObjectGetComm((PetscObject)v,&comm);

145:   VecGetBlockSize(v,&bs);
146:   if (start < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative start %D",start);
147:   else if (start >= bs) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Start of stride subvector (%D) is too large for stride\n Have you set the vector blocksize (%D) correctly with VecSetBlockSize()?",start,bs);
148:   x += start;

150:   if (ntype == NORM_2) {
151:     PetscScalar sum = 0.0;
152:     for (i=0; i<n; i+=bs) sum += x[i]*(PetscConj(x[i]));
153:     tnorm = PetscRealPart(sum);
154:     MPIU_Allreduce(&tnorm,nrm,1,MPIU_REAL,MPIU_SUM,comm);
155:     *nrm  = PetscSqrtReal(*nrm);
156:   } else if (ntype == NORM_1) {
157:     tnorm = 0.0;
158:     for (i=0; i<n; i+=bs) tnorm += PetscAbsScalar(x[i]);
159:     MPIU_Allreduce(&tnorm,nrm,1,MPIU_REAL,MPIU_SUM,comm);
160:   } else if (ntype == NORM_INFINITY) {
161:     PetscReal tmp;
162:     tnorm = 0.0;

164:     for (i=0; i<n; i+=bs) {
165:       if ((tmp = PetscAbsScalar(x[i])) > tnorm) tnorm = tmp;
166:       /* check special case of tmp == NaN */
167:       if (tmp != tmp) {tnorm = tmp; break;}
168:     }
169:     MPIU_Allreduce(&tnorm,nrm,1,MPIU_REAL,MPIU_MAX,comm);
170:   } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown norm type");
171:   VecRestoreArrayRead(v,&x);
172:   return(0);
173: }

175: /*@
176:    VecStrideMax - Computes the maximum of subvector of a vector defined
177:    by a starting point and a stride and optionally its location.

179:    Collective on Vec

181:    Input Parameter:
182: +  v - the vector
183: -  start - starting point of the subvector (defined by a stride)

185:    Output Parameter:
186: +  index - the location where the maximum occurred  (pass NULL if not required)
187: -  nrm - the maximum value in the subvector

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

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

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

200:    Level: advanced


203: .seealso: VecMax(), VecStrideNorm(), VecStrideGather(), VecStrideScatter(), VecStrideMin()
204: @*/
205: PetscErrorCode  VecStrideMax(Vec v,PetscInt start,PetscInt *idex,PetscReal *nrm)
206: {
207:   PetscErrorCode    ierr;
208:   PetscInt          i,n,bs,id;
209:   const PetscScalar *x;
210:   PetscReal         max,tmp;
211:   MPI_Comm          comm;


217:   VecGetLocalSize(v,&n);
218:   VecGetArrayRead(v,&x);
219:   PetscObjectGetComm((PetscObject)v,&comm);

221:   VecGetBlockSize(v,&bs);
222:   if (start < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative start %D",start);
223:   else if (start >= bs) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Start of stride subvector (%D) is too large for stride\n Have you set the vector blocksize (%D) correctly with VecSetBlockSize()?",start,bs);
224:   x += start;

226:   id = -1;
227:   if (!n) max = PETSC_MIN_REAL;
228:   else {
229:     id  = 0;
230:     max = PetscRealPart(x[0]);
231:     for (i=bs; i<n; i+=bs) {
232:       if ((tmp = PetscRealPart(x[i])) > max) { max = tmp; id = i;}
233:     }
234:   }
235:   VecRestoreArrayRead(v,&x);

237:   if (!idex) {
238:     MPIU_Allreduce(&max,nrm,1,MPIU_REAL,MPIU_MAX,comm);
239:   } else {
240:     PetscReal in[2],out[2];
241:     PetscInt  rstart;

243:     VecGetOwnershipRange(v,&rstart,NULL);
244:     in[0] = max;
245:     in[1] = rstart+id+start;
246:     MPIU_Allreduce(in,out,2,MPIU_REAL,MPIU_MAXINDEX_OP,PetscObjectComm((PetscObject)v));
247:     *nrm  = out[0];
248:     *idex = (PetscInt)out[1];
249:   }
250:   return(0);
251: }

253: /*@
254:    VecStrideMin - Computes the minimum of subvector of a vector defined
255:    by a starting point and a stride and optionally its location.

257:    Collective on Vec

259:    Input Parameter:
260: +  v - the vector
261: -  start - starting point of the subvector (defined by a stride)

263:    Output Parameter:
264: +  idex - the location where the minimum occurred. (pass NULL if not required)
265: -  nrm - the minimum value in the subvector

267:    Level: advanced

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

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

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


281: .seealso: VecMin(), VecStrideNorm(), VecStrideGather(), VecStrideScatter(), VecStrideMax()
282: @*/
283: PetscErrorCode  VecStrideMin(Vec v,PetscInt start,PetscInt *idex,PetscReal *nrm)
284: {
285:   PetscErrorCode    ierr;
286:   PetscInt          i,n,bs,id;
287:   const PetscScalar *x;
288:   PetscReal         min,tmp;
289:   MPI_Comm          comm;


295:   VecGetLocalSize(v,&n);
296:   VecGetArrayRead(v,&x);
297:   PetscObjectGetComm((PetscObject)v,&comm);

299:   VecGetBlockSize(v,&bs);
300:   if (start < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative start %D",start);
301:   else if (start >= bs) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Start of stride subvector (%D) is too large for stride\nHave you set the vector blocksize (%D) correctly with VecSetBlockSize()?",start,bs);
302:   x += start;

304:   id = -1;
305:   if (!n) min = PETSC_MAX_REAL;
306:   else {
307:     id = 0;
308:     min = PetscRealPart(x[0]);
309:     for (i=bs; i<n; i+=bs) {
310:       if ((tmp = PetscRealPart(x[i])) < min) { min = tmp; id = i;}
311:     }
312:   }
313:   VecRestoreArrayRead(v,&x);

315:   if (!idex) {
316:     MPIU_Allreduce(&min,nrm,1,MPIU_REAL,MPIU_MIN,comm);
317:   } else {
318:     PetscReal in[2],out[2];
319:     PetscInt  rstart;

321:     VecGetOwnershipRange(v,&rstart,NULL);
322:     in[0] = min;
323:     in[1] = rstart+id;
324:     MPIU_Allreduce(in,out,2,MPIU_REAL,MPIU_MININDEX_OP,PetscObjectComm((PetscObject)v));
325:     *nrm  = out[0];
326:     *idex = (PetscInt)out[1];
327:   }
328:   return(0);
329: }

331: /*@
332:    VecStrideScaleAll - Scales the subvectors of a vector defined
333:    by a starting point and a stride.

335:    Logically Collective on Vec

337:    Input Parameter:
338: +  v - the vector
339: -  scales - values to multiply each subvector entry by

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

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


348:    Level: advanced


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

362:   VecSetErrorIfLocked(v,1);

364:   VecGetLocalSize(v,&n);
365:   VecGetArray(v,&x);
366:   VecGetBlockSize(v,&bs);

368:   /* need to provide optimized code for each bs */
369:   for (i=0; i<n; i+=bs) {
370:     for (j=0; j<bs; j++) x[i+j] *= scales[j];
371:   }
372:   VecRestoreArray(v,&x);
373:   return(0);
374: }

376: /*@
377:    VecStrideNormAll - Computes the norms of subvectors of a vector defined
378:    by a starting point and a stride.

380:    Collective on Vec

382:    Input Parameter:
383: +  v - the vector
384: -  ntype - type of norm, one of NORM_1, NORM_2, NORM_INFINITY

386:    Output Parameter:
387: .  nrm - the norms

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

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

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

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

400:    Level: advanced


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

416:   VecGetLocalSize(v,&n);
417:   VecGetArrayRead(v,&x);
418:   PetscObjectGetComm((PetscObject)v,&comm);

420:   VecGetBlockSize(v,&bs);
421:   if (bs > 128) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Currently supports only blocksize up to 128");

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

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

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

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

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

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

462:    Collective on Vec

464:    Input Parameter:
465: .  v - the vector

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

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

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

478:    Level: advanced


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

494:   if (idex) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No support yet for returning index; send mail to petsc-maint@mcs.anl.gov asking for it");
495:   VecGetLocalSize(v,&n);
496:   VecGetArrayRead(v,&x);
497:   PetscObjectGetComm((PetscObject)v,&comm);

499:   VecGetBlockSize(v,&bs);
500:   if (bs > 128) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Currently supports only blocksize up to 128");

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

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

515:   VecRestoreArrayRead(v,&x);
516:   return(0);
517: }

519: /*@
520:    VecStrideMinAll - Computes the minimum of subvector of a vector defined
521:    by a starting point and a stride and optionally its location.

523:    Collective on Vec

525:    Input Parameter:
526: .  v - the vector

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

533:    Level: advanced

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

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


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

555:   if (idex) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No support yet for returning index; send mail to petsc-maint@mcs.anl.gov asking for it");
556:   VecGetLocalSize(v,&n);
557:   VecGetArrayRead(v,&x);
558:   PetscObjectGetComm((PetscObject)v,&comm);

560:   VecGetBlockSize(v,&bs);
561:   if (bs > 128) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Currently supports only blocksize up to 128");

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

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

576:   VecRestoreArrayRead(v,&x);
577:   return(0);
578: }

580: /*----------------------------------------------------------------------------------------------*/
581: /*@
582:    VecStrideGatherAll - Gathers all the single components from a multi-component vector into
583:    separate vectors.

585:    Collective on Vec

587:    Input Parameter:
588: +  v - the vector
589: -  addv - one of ADD_VALUES,INSERT_VALUES,MAX_VALUES

591:    Output Parameter:
592: .  s - the location where the subvectors are stored

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

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

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

605:    Not optimized; could be easily

607:    Level: advanced

609: .seealso: VecStrideNorm(), VecStrideScatter(), VecStrideMin(), VecStrideMax(), VecStrideGather(),
610:           VecStrideScatterAll()
611: @*/
612: PetscErrorCode  VecStrideGatherAll(Vec v,Vec s[],InsertMode addv)
613: {
614:   PetscErrorCode    ierr;
615:   PetscInt          i,n,n2,bs,j,k,*bss = NULL,nv,jj,nvc;
616:   PetscScalar       **y;
617:   const PetscScalar *x;

623:   VecGetLocalSize(v,&n);
624:   VecGetLocalSize(s[0],&n2);
625:   VecGetArrayRead(v,&x);
626:   VecGetBlockSize(v,&bs);
627:   if (bs <= 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Input vector does not have a valid blocksize set");

629:   PetscMalloc2(bs,&y,bs,&bss);
630:   nv   = 0;
631:   nvc  = 0;
632:   for (i=0; i<bs; i++) {
633:     VecGetBlockSize(s[i],&bss[i]);
634:     if (bss[i] < 1) bss[i] = 1; /* if user never set it then assume 1  Re: [PETSC #8241] VecStrideGatherAll */
635:     VecGetArray(s[i],&y[i]);
636:     nvc += bss[i];
637:     nv++;
638:     if (nvc > bs)  SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Number of subvectors in subvectors > number of vectors in main vector");
639:     if (nvc == bs) break;
640:   }

642:   n =  n/bs;

644:   jj = 0;
645:   if (addv == INSERT_VALUES) {
646:     for (j=0; j<nv; j++) {
647:       for (k=0; k<bss[j]; k++) {
648:         for (i=0; i<n; i++) y[j][i*bss[j] + k] = x[bs*i+jj+k];
649:       }
650:       jj += bss[j];
651:     }
652:   } else if (addv == ADD_VALUES) {
653:     for (j=0; j<nv; j++) {
654:       for (k=0; k<bss[j]; k++) {
655:         for (i=0; i<n; i++) y[j][i*bss[j] + k] += x[bs*i+jj+k];
656:       }
657:       jj += bss[j];
658:     }
659: #if !defined(PETSC_USE_COMPLEX)
660:   } else if (addv == MAX_VALUES) {
661:     for (j=0; j<nv; j++) {
662:       for (k=0; k<bss[j]; k++) {
663:         for (i=0; i<n; i++) y[j][i*bss[j] + k] = PetscMax(y[j][i*bss[j] + k],x[bs*i+jj+k]);
664:       }
665:       jj += bss[j];
666:     }
667: #endif
668:   } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown insert type");

670:   VecRestoreArrayRead(v,&x);
671:   for (i=0; i<nv; i++) {
672:     VecRestoreArray(s[i],&y[i]);
673:   }

675:   PetscFree2(y,bss);
676:   return(0);
677: }

679: /*@
680:    VecStrideScatterAll - Scatters all the single components from separate vectors into
681:      a multi-component vector.

683:    Collective on Vec

685:    Input Parameter:
686: +  s - the location where the subvectors are stored
687: -  addv - one of ADD_VALUES,INSERT_VALUES,MAX_VALUES

689:    Output Parameter:
690: .  v - the multicomponent vector

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

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

699:    Not optimized; could be easily

701:    Level: advanced

703: .seealso: VecStrideNorm(), VecStrideScatter(), VecStrideMin(), VecStrideMax(), VecStrideGather(),
704:           VecStrideScatterAll()
705: @*/
706: PetscErrorCode  VecStrideScatterAll(Vec s[],Vec v,InsertMode addv)
707: {
708:   PetscErrorCode    ierr;
709:   PetscInt          i,n,n2,bs,j,jj,k,*bss = NULL,nv,nvc;
710:   PetscScalar       *x;
711:   PetscScalar const **y;

717:   VecGetLocalSize(v,&n);
718:   VecGetLocalSize(s[0],&n2);
719:   VecGetArray(v,&x);
720:   VecGetBlockSize(v,&bs);
721:   if (bs <= 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Input vector does not have a valid blocksize set");

723:   PetscMalloc2(bs,(PetscScalar***)&y,bs,&bss);
724:   nv   = 0;
725:   nvc  = 0;
726:   for (i=0; i<bs; i++) {
727:     VecGetBlockSize(s[i],&bss[i]);
728:     if (bss[i] < 1) bss[i] = 1; /* if user never set it then assume 1  Re: [PETSC #8241] VecStrideGatherAll */
729:     VecGetArrayRead(s[i],&y[i]);
730:     nvc += bss[i];
731:     nv++;
732:     if (nvc > bs) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Number of subvectors in subvectors > number of vectors in main vector");
733:     if (nvc == bs) break;
734:   }

736:   n =  n/bs;

738:   jj = 0;
739:   if (addv == INSERT_VALUES) {
740:     for (j=0; j<nv; j++) {
741:       for (k=0; k<bss[j]; k++) {
742:         for (i=0; i<n; i++) x[bs*i+jj+k] = y[j][i*bss[j] + k];
743:       }
744:       jj += bss[j];
745:     }
746:   } else if (addv == ADD_VALUES) {
747:     for (j=0; j<nv; j++) {
748:       for (k=0; k<bss[j]; k++) {
749:         for (i=0; i<n; i++) x[bs*i+jj+k] += y[j][i*bss[j] + k];
750:       }
751:       jj += bss[j];
752:     }
753: #if !defined(PETSC_USE_COMPLEX)
754:   } else if (addv == MAX_VALUES) {
755:     for (j=0; j<nv; j++) {
756:       for (k=0; k<bss[j]; k++) {
757:         for (i=0; i<n; i++) x[bs*i+jj+k] = PetscMax(x[bs*i+jj+k],y[j][i*bss[j] + k]);
758:       }
759:       jj += bss[j];
760:     }
761: #endif
762:   } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown insert type");

764:   VecRestoreArray(v,&x);
765:   for (i=0; i<nv; i++) {
766:     VecRestoreArrayRead(s[i],&y[i]);
767:   }
768:   PetscFree2(*(PetscScalar***)&y,bss);
769:   return(0);
770: }

772: /*@
773:    VecStrideGather - Gathers a single component from a multi-component vector into
774:    another vector.

776:    Collective on Vec

778:    Input Parameter:
779: +  v - the vector
780: .  start - starting point of the subvector (defined by a stride)
781: -  addv - one of ADD_VALUES,INSERT_VALUES,MAX_VALUES

783:    Output Parameter:
784: .  s - the location where the subvector is stored

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

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

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

796:    Not optimized; could be easily

798:    Level: advanced

800: .seealso: VecStrideNorm(), VecStrideScatter(), VecStrideMin(), VecStrideMax(), VecStrideGatherAll(),
801:           VecStrideScatterAll()
802: @*/
803: PetscErrorCode  VecStrideGather(Vec v,PetscInt start,Vec s,InsertMode addv)
804: {

810:   if (start < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative start %D",start);
811:   if (start >= v->map->bs) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Start of stride subvector (%D) is too large for stride\n Have you set the vector blocksize (%D) correctly with VecSetBlockSize()?",start,v->map->bs);
812:   if (!v->ops->stridegather) SETERRQ(PetscObjectComm((PetscObject)s),PETSC_ERR_SUP,"Not implemented for this Vec class");
813:   (*v->ops->stridegather)(v,start,s,addv);
814:   return(0);
815: }

817: /*@
818:    VecStrideScatter - Scatters a single component from a vector into a multi-component vector.

820:    Collective on Vec

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

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

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

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

837:    Not optimized; could be easily

839:    Level: advanced

841: .seealso: VecStrideNorm(), VecStrideGather(), VecStrideMin(), VecStrideMax(), VecStrideGatherAll(),
842:           VecStrideScatterAll(), VecStrideSubSetScatter(), VecStrideSubSetGather()
843: @*/
844: PetscErrorCode  VecStrideScatter(Vec s,PetscInt start,Vec v,InsertMode addv)
845: {

851:   if (start < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative start %D",start);
852:   if (start >= v->map->bs) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Start of stride subvector (%D) is too large for stride\n Have you set the vector blocksize (%D) correctly with VecSetBlockSize()?",start,v->map->bs);
853:   if (!v->ops->stridescatter) SETERRQ(PetscObjectComm((PetscObject)s),PETSC_ERR_SUP,"Not implemented for this Vec class");
854:   (*v->ops->stridescatter)(s,start,v,addv);
855:   return(0);
856: }

858: /*@
859:    VecStrideSubSetGather - Gathers a subset of components from a multi-component vector into
860:    another vector.

862:    Collective on Vec

864:    Input Parameter:
865: +  v - the vector
866: .  nidx - the number of indices
867: .  idxv - the indices of the components 0 <= idxv[0] ...idxv[nidx-1] < bs(v), they need not be sorted
868: .  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
869: -  addv - one of ADD_VALUES,INSERT_VALUES,MAX_VALUES

871:    Output Parameter:
872: .  s - the location where the subvector is stored

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


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

881:    Not optimized; could be easily

883:    Level: advanced

885: .seealso: VecStrideNorm(), VecStrideScatter(), VecStrideGather(), VecStrideSubSetScatter(), VecStrideMin(), VecStrideMax(), VecStrideGatherAll(),
886:           VecStrideScatterAll()
887: @*/
888: PetscErrorCode  VecStrideSubSetGather(Vec v,PetscInt nidx,const PetscInt idxv[],const PetscInt idxs[],Vec s,InsertMode addv)
889: {

895:   if (nidx == PETSC_DETERMINE) nidx = s->map->bs;
896:   if (!v->ops->stridesubsetgather) SETERRQ(PetscObjectComm((PetscObject)s),PETSC_ERR_SUP,"Not implemented for this Vec class");
897:   (*v->ops->stridesubsetgather)(v,nidx,idxv,idxs,s,addv);
898:   return(0);
899: }

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

904:    Collective on Vec

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

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

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

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

922:    Not optimized; could be easily

924:    Level: advanced

926: .seealso: VecStrideNorm(), VecStrideGather(), VecStrideGather(), VecStrideSubSetGather(), VecStrideMin(), VecStrideMax(), VecStrideGatherAll(),
927:           VecStrideScatterAll()
928: @*/
929: PetscErrorCode  VecStrideSubSetScatter(Vec s,PetscInt nidx,const PetscInt idxs[],const PetscInt idxv[],Vec v,InsertMode addv)
930: {

936:   if (nidx == PETSC_DETERMINE) nidx = s->map->bs;
937:   if (!v->ops->stridesubsetscatter) SETERRQ(PetscObjectComm((PetscObject)s),PETSC_ERR_SUP,"Not implemented for this Vec class");
938:   (*v->ops->stridesubsetscatter)(s,nidx,idxs,idxv,v,addv);
939:   return(0);
940: }

942: PetscErrorCode  VecStrideGather_Default(Vec v,PetscInt start,Vec s,InsertMode addv)
943: {
945:   PetscInt       i,n,bs,ns;
946:   const PetscScalar *x;
947:   PetscScalar       *y;

950:   VecGetLocalSize(v,&n);
951:   VecGetLocalSize(s,&ns);
952:   VecGetArrayRead(v,&x);
953:   VecGetArray(s,&y);

955:   bs = v->map->bs;
956:   if (n != ns*bs) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Subvector length * blocksize %D not correct for gather from original vector %D",ns*bs,n);
957:   x += start;
958:   n  =  n/bs;

960:   if (addv == INSERT_VALUES) {
961:     for (i=0; i<n; i++) y[i] = x[bs*i];
962:   } else if (addv == ADD_VALUES) {
963:     for (i=0; i<n; i++) y[i] += x[bs*i];
964: #if !defined(PETSC_USE_COMPLEX)
965:   } else if (addv == MAX_VALUES) {
966:     for (i=0; i<n; i++) y[i] = PetscMax(y[i],x[bs*i]);
967: #endif
968:   } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown insert type");

970:   VecRestoreArrayRead(v,&x);
971:   VecRestoreArray(s,&y);
972:   return(0);
973: }

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

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

988:   bs = v->map->bs;
989:   if (n != ns*bs) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Subvector length * blocksize %D not correct for scatter to multicomponent vector %D",ns*bs,n);
990:   x += start;
991:   n  =  n/bs;

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

1003:   VecRestoreArray(v,&x);
1004:   VecRestoreArrayRead(s,&y);
1005:   return(0);
1006: }

1008: PetscErrorCode  VecStrideSubSetGather_Default(Vec v,PetscInt nidx,const PetscInt idxv[],const PetscInt idxs[],Vec s,InsertMode addv)
1009: {
1010:   PetscErrorCode    ierr;
1011:   PetscInt          i,j,n,bs,bss,ns;
1012:   const PetscScalar *x;
1013:   PetscScalar       *y;

1016:   VecGetLocalSize(v,&n);
1017:   VecGetLocalSize(s,&ns);
1018:   VecGetArrayRead(v,&x);
1019:   VecGetArray(s,&y);

1021:   bs  = v->map->bs;
1022:   bss = s->map->bs;
1023:   n  =  n/bs;

1025:   if (PetscDefined(USE_DEBUG)) {
1026:     if (n != ns/bss) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Incompatible layout of vectors");
1027:     for (j=0; j<nidx; j++) {
1028:       if (idxv[j] < 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"idx[%D] %D is negative",j,idxv[j]);
1029:       if (idxv[j] >= bs) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"idx[%D] %D is greater than or equal to vector blocksize %D",j,idxv[j],bs);
1030:     }
1031:     if (!idxs && bss != nidx) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Must provide idxs when not gathering into all locations");
1032:   }

1034:   if (addv == INSERT_VALUES) {
1035:     if (!idxs) {
1036:       for (i=0; i<n; i++) {
1037:         for (j=0; j<bss; j++) y[bss*i+j] = x[bs*i+idxv[j]];
1038:       }
1039:     } else {
1040:       for (i=0; i<n; i++) {
1041:         for (j=0; j<bss; j++) y[bss*i+idxs[j]] = x[bs*i+idxv[j]];
1042:       }
1043:     }
1044:   } else if (addv == ADD_VALUES) {
1045:     if (!idxs) {
1046:       for (i=0; i<n; i++) {
1047:         for (j=0; j<bss; j++) y[bss*i+j] += x[bs*i+idxv[j]];
1048:       }
1049:     } else {
1050:       for (i=0; i<n; i++) {
1051:         for (j=0; j<bss; j++) y[bss*i+idxs[j]] += x[bs*i+idxv[j]];
1052:       }
1053:     }
1054: #if !defined(PETSC_USE_COMPLEX)
1055:   } else if (addv == MAX_VALUES) {
1056:     if (!idxs) {
1057:       for (i=0; i<n; i++) {
1058:         for (j=0; j<bss; j++) y[bss*i+j] = PetscMax(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]] = PetscMax(y[bss*i+idxs[j]],x[bs*i+idxv[j]]);
1063:       }
1064:     }
1065: #endif
1066:   } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown insert type");

1068:   VecRestoreArrayRead(v,&x);
1069:   VecRestoreArray(s,&y);
1070:   return(0);
1071: }

1073: PetscErrorCode  VecStrideSubSetScatter_Default(Vec s,PetscInt nidx,const PetscInt idxs[],const PetscInt idxv[],Vec v,InsertMode addv)
1074: {
1075:   PetscErrorCode    ierr;
1076:   PetscInt          j,i,n,bs,ns,bss;
1077:   PetscScalar       *x;
1078:   const PetscScalar *y;

1081:   VecGetLocalSize(v,&n);
1082:   VecGetLocalSize(s,&ns);
1083:   VecGetArray(v,&x);
1084:   VecGetArrayRead(s,&y);

1086:   bs  = v->map->bs;
1087:   bss = s->map->bs;
1088:   n  =  n/bs;

1090:   if (PetscDefined(USE_DEBUG)) {
1091:     if (n != ns/bss) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Incompatible layout of vectors");
1092:     for (j=0; j<bss; j++) {
1093:       if (idxs) {
1094:         if (idxs[j] < 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"idx[%D] %D is negative",j,idxs[j]);
1095:         if (idxs[j] >= bs) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"idx[%D] %D is greater than or equal to vector blocksize %D",j,idxs[j],bs);
1096:       }
1097:     }
1098:     if (!idxs && bss != nidx) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Must provide idxs when not scattering from all locations");
1099:   }

1101:   if (addv == INSERT_VALUES) {
1102:     if (!idxs) {
1103:       for (i=0; i<n; i++) {
1104:         for (j=0; j<bss; j++) x[bs*i + idxv[j]] = y[bss*i+j];
1105:       }
1106:     } else {
1107:       for (i=0; i<n; i++) {
1108:         for (j=0; j<bss; j++) x[bs*i + idxv[j]] = y[bss*i+idxs[j]];
1109:       }
1110:     }
1111:   } else if (addv == ADD_VALUES) {
1112:     if (!idxs) {
1113:       for (i=0; i<n; i++) {
1114:         for (j=0; j<bss; j++) x[bs*i + idxv[j]] += y[bss*i+j];
1115:       }
1116:     } else {
1117:       for (i=0; i<n; i++) {
1118:         for (j=0; j<bss; j++) x[bs*i + idxv[j]] += y[bss*i+idxs[j]];
1119:       }
1120:     }
1121: #if !defined(PETSC_USE_COMPLEX)
1122:   } else if (addv == MAX_VALUES) {
1123:     if (!idxs) {
1124:       for (i=0; i<n; i++) {
1125:         for (j=0; j<bss; j++) x[bs*i + idxv[j]] = PetscMax(y[bss*i+j],x[bs*i + idxv[j]]);
1126:       }
1127:     } else {
1128:       for (i=0; i<n; i++) {
1129:         for (j=0; j<bss; j++) x[bs*i + idxv[j]] = PetscMax(y[bss*i+idxs[j]],x[bs*i + idxv[j]]);
1130:       }
1131:     }
1132: #endif
1133:   } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown insert type");

1135:   VecRestoreArray(v,&x);
1136:   VecRestoreArrayRead(s,&y);
1137:   return(0);
1138: }

1140: PetscErrorCode VecReciprocal_Default(Vec v)
1141: {
1143:   PetscInt       i,n;
1144:   PetscScalar    *x;

1147:   VecGetLocalSize(v,&n);
1148:   VecGetArray(v,&x);
1149:   for (i=0; i<n; i++) {
1150:     if (x[i] != (PetscScalar)0.0) x[i] = (PetscScalar)1.0/x[i];
1151:   }
1152:   VecRestoreArray(v,&x);
1153:   return(0);
1154: }

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

1159:   Not collective

1161:   Input Parameter:
1162: . v - The vector

1164:   Output Parameter:
1165: . v - The vector of exponents

1167:   Level: beginner

1169: .seealso:  VecLog(), VecAbs(), VecSqrtAbs(), VecReciprocal()

1171: @*/
1172: PetscErrorCode  VecExp(Vec v)
1173: {
1174:   PetscScalar    *x;
1175:   PetscInt       i, n;

1180:   if (v->ops->exp) {
1181:     (*v->ops->exp)(v);
1182:   } else {
1183:     VecGetLocalSize(v, &n);
1184:     VecGetArray(v, &x);
1185:     for (i = 0; i < n; i++) x[i] = PetscExpScalar(x[i]);
1186:     VecRestoreArray(v, &x);
1187:   }
1188:   return(0);
1189: }

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

1194:   Not collective

1196:   Input Parameter:
1197: . v - The vector

1199:   Output Parameter:
1200: . v - The vector of logs

1202:   Level: beginner

1204: .seealso:  VecExp(), VecAbs(), VecSqrtAbs(), VecReciprocal()

1206: @*/
1207: PetscErrorCode  VecLog(Vec v)
1208: {
1209:   PetscScalar    *x;
1210:   PetscInt       i, n;

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

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

1229:   Not collective

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

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

1237:   Level: beginner

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

1241: .seealso: VecLog(), VecExp(), VecReciprocal(), VecAbs()

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

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

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

1266:   Collective on Vec

1268:   Input Parameter:
1269: + s - first vector
1270: - t - second vector

1272:   Output Parameter:
1273: + dp - s'conj(t)
1274: - nm - t'conj(t)

1276:   Level: advanced

1278:   Notes:
1279:     conj(x) is the complex conjugate of x when x is complex


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

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

1300:   if (s->map->N != t->map->N) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Incompatible vector global lengths");
1301:   if (s->map->n != t->map->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Incompatible vector local lengths");

1303:   PetscLogEventBegin(VEC_DotNorm2,s,t,0,0);
1304:   if (s->ops->dotnorm2) {
1305:     (*s->ops->dotnorm2)(s,t,dp,&dpx);
1306:     *nm  = PetscRealPart(dpx);
1307:   } else {
1308:     VecGetLocalSize(s, &n);
1309:     VecGetArrayRead(s, &sx);
1310:     VecGetArrayRead(t, &tx);

1312:     for (i = 0; i<n; i++) {
1313:       dpx += sx[i]*PetscConj(tx[i]);
1314:       nmx += tx[i]*PetscConj(tx[i]);
1315:     }
1316:     work[0] = dpx;
1317:     work[1] = nmx;

1319:     MPIU_Allreduce(work,sum,2,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)s));
1320:     *dp  = sum[0];
1321:     *nm  = PetscRealPart(sum[1]);

1323:     VecRestoreArrayRead(t, &tx);
1324:     VecRestoreArrayRead(s, &sx);
1325:     PetscLogFlops(4.0*n);
1326:   }
1327:   PetscLogEventEnd(VEC_DotNorm2,s,t,0,0);
1328:   return(0);
1329: }

1331: /*@
1332:    VecSum - Computes the sum of all the components of a vector.

1334:    Collective on Vec

1336:    Input Parameter:
1337: .  v - the vector

1339:    Output Parameter:
1340: .  sum - the result

1342:    Level: beginner

1344: .seealso: VecNorm()
1345: @*/
1346: PetscErrorCode  VecSum(Vec v,PetscScalar *sum)
1347: {
1348:   PetscErrorCode    ierr;
1349:   PetscInt          i,n;
1350:   const PetscScalar *x;
1351:   PetscScalar       lsum = 0.0;

1356:   VecGetLocalSize(v,&n);
1357:   VecGetArrayRead(v,&x);
1358:   for (i=0; i<n; i++) lsum += x[i];
1359:   MPIU_Allreduce(&lsum,sum,1,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)v));
1360:   VecRestoreArrayRead(v,&x);
1361:   return(0);
1362: }

1364: /*@
1365:    VecImaginaryPart - Replaces a complex vector with its imginary part

1367:    Collective on Vec

1369:    Input Parameter:
1370: .  v - the vector

1372:    Level: beginner

1374: .seealso: VecNorm(), VecRealPart()
1375: @*/
1376: PetscErrorCode  VecImaginaryPart(Vec v)
1377: {
1378:   PetscErrorCode    ierr;
1379:   PetscInt          i,n;
1380:   PetscScalar       *x;

1384:   VecGetLocalSize(v,&n);
1385:   VecGetArray(v,&x);
1386:   for (i=0; i<n; i++) x[i] = PetscImaginaryPart(x[i]);
1387:   VecRestoreArray(v,&x);
1388:   return(0);
1389: }

1391: /*@
1392:    VecRealPart - Replaces a complex vector with its real part

1394:    Collective on Vec

1396:    Input Parameter:
1397: .  v - the vector

1399:    Level: beginner

1401: .seealso: VecNorm(), VecImaginaryPart()
1402: @*/
1403: PetscErrorCode  VecRealPart(Vec v)
1404: {
1405:   PetscErrorCode    ierr;
1406:   PetscInt          i,n;
1407:   PetscScalar       *x;

1411:   VecGetLocalSize(v,&n);
1412:   VecGetArray(v,&x);
1413:   for (i=0; i<n; i++) x[i] = PetscRealPart(x[i]);
1414:   VecRestoreArray(v,&x);
1415:   return(0);
1416: }

1418: /*@
1419:    VecShift - Shifts all of the components of a vector by computing
1420:    x[i] = x[i] + shift.

1422:    Logically Collective on Vec

1424:    Input Parameters:
1425: +  v - the vector
1426: -  shift - the shift

1428:    Output Parameter:
1429: .  v - the shifted vector

1431:    Level: intermediate

1433: @*/
1434: PetscErrorCode  VecShift(Vec v,PetscScalar shift)
1435: {
1437:   PetscInt       i,n;
1438:   PetscScalar    *x;

1443:   VecSetErrorIfLocked(v,1);
1444:   if (shift == 0.0) return(0);

1446:   if (v->ops->shift) {
1447:     (*v->ops->shift)(v,shift);
1448:   } else {
1449:     VecGetLocalSize(v,&n);
1450:     VecGetArray(v,&x);
1451:     for (i=0; i<n; i++) x[i] += shift;
1452:     VecRestoreArray(v,&x);
1453:   }
1454:   return(0);
1455: }

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

1460:    Logically Collective on Vec

1462:    Input Parameters:
1463: .  v - the vector

1465:    Level: intermediate

1467: @*/
1468: PetscErrorCode  VecAbs(Vec v)
1469: {
1471:   PetscInt       i,n;
1472:   PetscScalar    *x;

1476:   VecSetErrorIfLocked(v,1);

1478:   if (v->ops->abs) {
1479:     (*v->ops->abs)(v);
1480:   } else {
1481:     VecGetLocalSize(v,&n);
1482:     VecGetArray(v,&x);
1483:     for (i=0; i<n; i++) x[i] = PetscAbsScalar(x[i]);
1484:     VecRestoreArray(v,&x);
1485:   }
1486:   return(0);
1487: }

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

1492:   Input Parameters:
1493: + vec   - The vector
1494: . order - The ordering
1495: - inv   - The flag for inverting the permutation

1497:   Level: beginner

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

1501: .seealso: MatPermute()
1502: @*/
1503: PetscErrorCode  VecPermute(Vec x, IS row, PetscBool inv)
1504: {
1505:   const PetscScalar *array;
1506:   PetscScalar       *newArray;
1507:   const PetscInt    *idx;
1508:   PetscInt          i,rstart,rend;
1509:   PetscErrorCode    ierr;

1512:   VecSetErrorIfLocked(x,1);

1514:   VecGetOwnershipRange(x,&rstart,&rend);
1515:   ISGetIndices(row, &idx);
1516:   VecGetArrayRead(x, &array);
1517:   PetscMalloc1(x->map->n, &newArray);
1518:   if (PetscDefined(USE_DEBUG)) {
1519:     for (i = 0; i < x->map->n; i++) {
1520:       if ((idx[i] < rstart) || (idx[i] >= rend)) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_CORRUPT, "Permutation index %D is out of bounds: %D", i, idx[i]);
1521:     }
1522:   }
1523:   if (!inv) {
1524:     for (i = 0; i < x->map->n; i++) newArray[i] = array[idx[i]-rstart];
1525:   } else {
1526:     for (i = 0; i < x->map->n; i++) newArray[idx[i]-rstart] = array[i];
1527:   }
1528:   VecRestoreArrayRead(x, &array);
1529:   ISRestoreIndices(row, &idx);
1530:   VecReplaceArray(x, newArray);
1531:   return(0);
1532: }

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

1539:    Collective on Vec

1541:    Input Parameters:
1542: +  vec1 - the first vector
1543: -  vec2 - the second vector

1545:    Output Parameter:
1546: .  flg - PETSC_TRUE if the vectors are equal; PETSC_FALSE otherwise.

1548:    Level: intermediate


1551: @*/
1552: PetscErrorCode  VecEqual(Vec vec1,Vec vec2,PetscBool  *flg)
1553: {
1554:   const PetscScalar  *v1,*v2;
1555:   PetscErrorCode     ierr;
1556:   PetscInt           n1,n2,N1,N2;
1557:   PetscBool          flg1;

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

1586: /*@
1587:    VecUniqueEntries - Compute the number of unique entries, and those entries

1589:    Collective on Vec

1591:    Input Parameter:
1592: .  vec - the vector

1594:    Output Parameters:
1595: +  n - The number of unique entries
1596: -  e - The entries

1598:    Level: intermediate

1600: @*/
1601: PetscErrorCode  VecUniqueEntries(Vec vec, PetscInt *n, PetscScalar **e)
1602: {
1603:   const PetscScalar *v;
1604:   PetscScalar       *tmp, *vals;
1605:   PetscMPIInt       *N, *displs, l;
1606:   PetscInt          ng, m, i, j, p;
1607:   PetscMPIInt       size;
1608:   PetscErrorCode    ierr;

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