Actual source code: vinv.c

petsc-3.6.0 2015-06-09
Report Typos and Errors
  2: /*
  3:      Some useful vector utility functions.
  4: */
  5: #include <../src/vec/vec/impls/mpi/pvecimpl.h>          /*I "petscvec.h" I*/
  6: extern MPI_Op VecMax_Local_Op;
  7: extern MPI_Op VecMin_Local_Op;

 11: /*@
 12:    VecStrideSet - Sets a subvector of a vector defined
 13:    by a starting point and a stride with a given value

 15:    Logically Collective on Vec

 17:    Input Parameter:
 18: +  v - the vector
 19: .  start - starting point of the subvector (defined by a stride)
 20: -  s - value to set for each entry in that subvector

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

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

 28:    Level: advanced

 30:    Concepts: scale^on stride of vector
 31:    Concepts: stride^scale

 33: .seealso: VecNorm(), VecStrideGather(), VecStrideScatter(), VecStrideMin(), VecStrideMax(), VecStrideScale()
 34: @*/
 35: PetscErrorCode  VecStrideSet(Vec v,PetscInt start,PetscScalar s)
 36: {
 38:   PetscInt       i,n,bs;
 39:   PetscScalar    *x;

 45:   VecLocked(v,1);
 46: 
 47:   VecGetLocalSize(v,&n);
 48:   VecGetArray(v,&x);
 49:   VecGetBlockSize(v,&bs);
 50:   if (start < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative start %D",start);
 51:   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);
 52:   x += start;

 54:   for (i=0; i<n; i+=bs) x[i] = s;
 55:   x -= start;

 57:   VecRestoreArray(v,&x);
 58:   return(0);
 59: }

 63: /*@
 64:    VecStrideScale - Scales a subvector of a vector defined
 65:    by a starting point and a stride.

 67:    Logically Collective on Vec

 69:    Input Parameter:
 70: +  v - the vector
 71: .  start - starting point of the subvector (defined by a stride)
 72: -  scale - value to multiply each subvector entry by

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

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

 80:    Level: advanced

 82:    Concepts: scale^on stride of vector
 83:    Concepts: stride^scale

 85: .seealso: VecNorm(), VecStrideGather(), VecStrideScatter(), VecStrideMin(), VecStrideMax(), VecStrideScale()
 86: @*/
 87: PetscErrorCode  VecStrideScale(Vec v,PetscInt start,PetscScalar scale)
 88: {
 90:   PetscInt       i,n,bs;
 91:   PetscScalar    *x;

 97:   VecLocked(v,1);
 98: 
 99:   VecGetLocalSize(v,&n);
100:   VecGetArray(v,&x);
101:   VecGetBlockSize(v,&bs);
102:   if (start < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative start %D",start);
103:   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);
104:   x += start;

106:   for (i=0; i<n; i+=bs) x[i] *= scale;
107:   x -= start;

109:   VecRestoreArray(v,&x);
110:   return(0);
111: }

115: /*@
116:    VecStrideNorm - Computes the norm of subvector of a vector defined
117:    by a starting point and a stride.

119:    Collective on Vec

121:    Input Parameter:
122: +  v - the vector
123: .  start - starting point of the subvector (defined by a stride)
124: -  ntype - type of norm, one of NORM_1, NORM_2, NORM_INFINITY

126:    Output Parameter:
127: .  norm - the norm

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

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

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

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

141:    Level: advanced

143:    Concepts: norm^on stride of vector
144:    Concepts: stride^norm

146: .seealso: VecNorm(), VecStrideGather(), VecStrideScatter(), VecStrideMin(), VecStrideMax()
147: @*/
148: PetscErrorCode  VecStrideNorm(Vec v,PetscInt start,NormType ntype,PetscReal *nrm)
149: {
150:   PetscErrorCode    ierr;
151:   PetscInt          i,n,bs;
152:   const PetscScalar *x;
153:   PetscReal         tnorm;
154:   MPI_Comm          comm;

159:   VecGetLocalSize(v,&n);
160:   VecGetArrayRead(v,&x);
161:   PetscObjectGetComm((PetscObject)v,&comm);

163:   VecGetBlockSize(v,&bs);
164:   if (start < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative start %D",start);
165:   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);
166:   x += start;

168:   if (ntype == NORM_2) {
169:     PetscScalar sum = 0.0;
170:     for (i=0; i<n; i+=bs) sum += x[i]*(PetscConj(x[i]));
171:     tnorm = PetscRealPart(sum);
172:     MPI_Allreduce(&tnorm,nrm,1,MPIU_REAL,MPIU_SUM,comm);
173:     *nrm  = PetscSqrtReal(*nrm);
174:   } else if (ntype == NORM_1) {
175:     tnorm = 0.0;
176:     for (i=0; i<n; i+=bs) tnorm += PetscAbsScalar(x[i]);
177:     MPI_Allreduce(&tnorm,nrm,1,MPIU_REAL,MPIU_SUM,comm);
178:   } else if (ntype == NORM_INFINITY) {
179:     PetscReal tmp;
180:     tnorm = 0.0;

182:     for (i=0; i<n; i+=bs) {
183:       if ((tmp = PetscAbsScalar(x[i])) > tnorm) tnorm = tmp;
184:       /* check special case of tmp == NaN */
185:       if (tmp != tmp) {tnorm = tmp; break;}
186:     }
187:     MPI_Allreduce(&tnorm,nrm,1,MPIU_REAL,MPIU_MAX,comm);
188:   } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown norm type");
189:   VecRestoreArrayRead(v,&x);
190:   return(0);
191: }

195: /*@
196:    VecStrideMax - Computes the maximum of subvector of a vector defined
197:    by a starting point and a stride and optionally its location.

199:    Collective on Vec

201:    Input Parameter:
202: +  v - the vector
203: -  start - starting point of the subvector (defined by a stride)

205:    Output Parameter:
206: +  index - the location where the maximum occurred  (pass NULL if not required)
207: -  nrm - the maximum value in the subvector

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

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

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

220:    Level: advanced

222:    Concepts: maximum^on stride of vector
223:    Concepts: stride^maximum

225: .seealso: VecMax(), VecStrideNorm(), VecStrideGather(), VecStrideScatter(), VecStrideMin()
226: @*/
227: PetscErrorCode  VecStrideMax(Vec v,PetscInt start,PetscInt *idex,PetscReal *nrm)
228: {
229:   PetscErrorCode    ierr;
230:   PetscInt          i,n,bs,id;
231:   const PetscScalar *x;
232:   PetscReal         max,tmp;
233:   MPI_Comm          comm;


239:   VecGetLocalSize(v,&n);
240:   VecGetArrayRead(v,&x);
241:   PetscObjectGetComm((PetscObject)v,&comm);

243:   VecGetBlockSize(v,&bs);
244:   if (start < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative start %D",start);
245:   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);
246:   x += start;

248:   id = -1;
249:   if (!n) max = PETSC_MIN_REAL;
250:   else {
251:     id  = 0;
252:     max = PetscRealPart(x[0]);
253:     for (i=bs; i<n; i+=bs) {
254:       if ((tmp = PetscRealPart(x[i])) > max) { max = tmp; id = i;}
255:     }
256:   }
257:   VecRestoreArrayRead(v,&x);

259:   if (!idex) {
260:     MPI_Allreduce(&max,nrm,1,MPIU_REAL,MPIU_MAX,comm);
261:   } else {
262:     PetscReal in[2],out[2];
263:     PetscInt  rstart;

265:     VecGetOwnershipRange(v,&rstart,NULL);
266:     in[0] = max;
267:     in[1] = rstart+id+start;
268:     MPI_Allreduce(in,out,2,MPIU_REAL,VecMax_Local_Op,PetscObjectComm((PetscObject)v));
269:     *nrm  = out[0];
270:     *idex = (PetscInt)out[1];
271:   }
272:   return(0);
273: }

277: /*@
278:    VecStrideMin - Computes the minimum of subvector of a vector defined
279:    by a starting point and a stride and optionally its location.

281:    Collective on Vec

283:    Input Parameter:
284: +  v - the vector
285: -  start - starting point of the subvector (defined by a stride)

287:    Output Parameter:
288: +  idex - the location where the minimum occurred. (pass NULL if not required)
289: -  nrm - the minimum value in the subvector

291:    Level: advanced

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

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

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

304:    Concepts: minimum^on stride of vector
305:    Concepts: stride^minimum

307: .seealso: VecMin(), VecStrideNorm(), VecStrideGather(), VecStrideScatter(), VecStrideMax()
308: @*/
309: PetscErrorCode  VecStrideMin(Vec v,PetscInt start,PetscInt *idex,PetscReal *nrm)
310: {
311:   PetscErrorCode    ierr;
312:   PetscInt          i,n,bs,id;
313:   const PetscScalar *x;
314:   PetscReal         min,tmp;
315:   MPI_Comm          comm;


321:   VecGetLocalSize(v,&n);
322:   VecGetArrayRead(v,&x);
323:   PetscObjectGetComm((PetscObject)v,&comm);

325:   VecGetBlockSize(v,&bs);
326:   if (start < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative start %D",start);
327:   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);
328:   x += start;

330:   id = -1;
331:   if (!n) min = PETSC_MAX_REAL;
332:   else {
333:     id = 0;
334:     min = PetscRealPart(x[0]);
335:     for (i=bs; i<n; i+=bs) {
336:       if ((tmp = PetscRealPart(x[i])) < min) { min = tmp; id = i;}
337:     }
338:   }
339:   VecRestoreArrayRead(v,&x);

341:   if (!idex) {
342:     MPI_Allreduce(&min,nrm,1,MPIU_REAL,MPIU_MIN,comm);
343:   } else {
344:     PetscReal in[2],out[2];
345:     PetscInt  rstart;

347:     VecGetOwnershipRange(v,&rstart,NULL);
348:     in[0] = min;
349:     in[1] = rstart+id;
350:     MPI_Allreduce(in,out,2,MPIU_REAL,VecMin_Local_Op,PetscObjectComm((PetscObject)v));
351:     *nrm  = out[0];
352:     *idex = (PetscInt)out[1];
353:   }
354:   return(0);
355: }

359: /*@
360:    VecStrideScaleAll - Scales the subvectors of a vector defined
361:    by a starting point and a stride.

363:    Logically Collective on Vec

365:    Input Parameter:
366: +  v - the vector
367: -  scales - values to multiply each subvector entry by

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

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


376:    Level: advanced

378:    Concepts: scale^on stride of vector
379:    Concepts: stride^scale

381: .seealso: VecNorm(), VecStrideScale(), VecScale(), VecStrideGather(), VecStrideScatter(), VecStrideMin(), VecStrideMax()
382: @*/
383: PetscErrorCode  VecStrideScaleAll(Vec v,const PetscScalar *scales)
384: {
386:   PetscInt       i,j,n,bs;
387:   PetscScalar    *x;

392:   VecLocked(v,1);
393: 
394:   VecGetLocalSize(v,&n);
395:   VecGetArray(v,&x);
396:   VecGetBlockSize(v,&bs);

398:   /* need to provide optimized code for each bs */
399:   for (i=0; i<n; i+=bs) {
400:     for (j=0; j<bs; j++) x[i+j] *= scales[j];
401:   }
402:   VecRestoreArray(v,&x);
403:   return(0);
404: }

408: /*@
409:    VecStrideNormAll - Computes the norms of subvectors of a vector defined
410:    by a starting point and a stride.

412:    Collective on Vec

414:    Input Parameter:
415: +  v - the vector
416: -  ntype - type of norm, one of NORM_1, NORM_2, NORM_INFINITY

418:    Output Parameter:
419: .  nrm - the norms

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

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

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

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

432:    Level: advanced

434:    Concepts: norm^on stride of vector
435:    Concepts: stride^norm

437: .seealso: VecNorm(), VecStrideGather(), VecStrideScatter(), VecStrideMin(), VecStrideMax()
438: @*/
439: PetscErrorCode  VecStrideNormAll(Vec v,NormType ntype,PetscReal nrm[])
440: {
441:   PetscErrorCode    ierr;
442:   PetscInt          i,j,n,bs;
443:   const PetscScalar *x;
444:   PetscReal         tnorm[128];
445:   MPI_Comm          comm;

450:   VecGetLocalSize(v,&n);
451:   VecGetArrayRead(v,&x);
452:   PetscObjectGetComm((PetscObject)v,&comm);

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

457:   if (ntype == NORM_2) {
458:     PetscScalar sum[128];
459:     for (j=0; j<bs; j++) sum[j] = 0.0;
460:     for (i=0; i<n; i+=bs) {
461:       for (j=0; j<bs; j++) sum[j] += x[i+j]*(PetscConj(x[i+j]));
462:     }
463:     for (j=0; j<bs; j++) tnorm[j]  = PetscRealPart(sum[j]);

465:     MPI_Allreduce(tnorm,nrm,bs,MPIU_REAL,MPIU_SUM,comm);
466:     for (j=0; j<bs; j++) nrm[j] = PetscSqrtReal(nrm[j]);
467:   } else if (ntype == NORM_1) {
468:     for (j=0; j<bs; j++) tnorm[j] = 0.0;

470:     for (i=0; i<n; i+=bs) {
471:       for (j=0; j<bs; j++) tnorm[j] += PetscAbsScalar(x[i+j]);
472:     }

474:     MPI_Allreduce(tnorm,nrm,bs,MPIU_REAL,MPIU_SUM,comm);
475:   } else if (ntype == NORM_INFINITY) {
476:     PetscReal tmp;
477:     for (j=0; j<bs; j++) tnorm[j] = 0.0;

479:     for (i=0; i<n; i+=bs) {
480:       for (j=0; j<bs; j++) {
481:         if ((tmp = PetscAbsScalar(x[i+j])) > tnorm[j]) tnorm[j] = tmp;
482:         /* check special case of tmp == NaN */
483:         if (tmp != tmp) {tnorm[j] = tmp; break;}
484:       }
485:     }
486:     MPI_Allreduce(tnorm,nrm,bs,MPIU_REAL,MPIU_MAX,comm);
487:   } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown norm type");
488:   VecRestoreArrayRead(v,&x);
489:   return(0);
490: }

494: /*@
495:    VecStrideMaxAll - Computes the maximums of subvectors of a vector defined
496:    by a starting point and a stride and optionally its location.

498:    Collective on Vec

500:    Input Parameter:
501: .  v - the vector

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

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

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

514:    Level: advanced

516:    Concepts: maximum^on stride of vector
517:    Concepts: stride^maximum

519: .seealso: VecMax(), VecStrideNorm(), VecStrideGather(), VecStrideScatter(), VecStrideMin()
520: @*/
521: PetscErrorCode  VecStrideMaxAll(Vec v,PetscInt idex[],PetscReal nrm[])
522: {
523:   PetscErrorCode    ierr;
524:   PetscInt          i,j,n,bs;
525:   const PetscScalar *x;
526:   PetscReal         max[128],tmp;
527:   MPI_Comm          comm;

532:   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");
533:   VecGetLocalSize(v,&n);
534:   VecGetArrayRead(v,&x);
535:   PetscObjectGetComm((PetscObject)v,&comm);

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

540:   if (!n) {
541:     for (j=0; j<bs; j++) max[j] = PETSC_MIN_REAL;
542:   } else {
543:     for (j=0; j<bs; j++) max[j] = PetscRealPart(x[j]);

545:     for (i=bs; i<n; i+=bs) {
546:       for (j=0; j<bs; j++) {
547:         if ((tmp = PetscRealPart(x[i+j])) > max[j]) max[j] = tmp;
548:       }
549:     }
550:   }
551:   MPI_Allreduce(max,nrm,bs,MPIU_REAL,MPIU_MAX,comm);

553:   VecRestoreArrayRead(v,&x);
554:   return(0);
555: }

559: /*@
560:    VecStrideMinAll - Computes the minimum of subvector of a vector defined
561:    by a starting point and a stride and optionally its location.

563:    Collective on Vec

565:    Input Parameter:
566: .  v - the vector

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

573:    Level: advanced

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

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

581:    Concepts: minimum^on stride of vector
582:    Concepts: stride^minimum

584: .seealso: VecMin(), VecStrideNorm(), VecStrideGather(), VecStrideScatter(), VecStrideMax()
585: @*/
586: PetscErrorCode  VecStrideMinAll(Vec v,PetscInt idex[],PetscReal nrm[])
587: {
588:   PetscErrorCode    ierr;
589:   PetscInt          i,n,bs,j;
590:   const PetscScalar *x;
591:   PetscReal         min[128],tmp;
592:   MPI_Comm          comm;

597:   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");
598:   VecGetLocalSize(v,&n);
599:   VecGetArrayRead(v,&x);
600:   PetscObjectGetComm((PetscObject)v,&comm);

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

605:   if (!n) {
606:     for (j=0; j<bs; j++) min[j] = PETSC_MAX_REAL;
607:   } else {
608:     for (j=0; j<bs; j++) min[j] = PetscRealPart(x[j]);

610:     for (i=bs; i<n; i+=bs) {
611:       for (j=0; j<bs; j++) {
612:         if ((tmp = PetscRealPart(x[i+j])) < min[j]) min[j] = tmp;
613:       }
614:     }
615:   }
616:   MPI_Allreduce(min,nrm,bs,MPIU_REAL,MPIU_MIN,comm);

618:   VecRestoreArrayRead(v,&x);
619:   return(0);
620: }

622: /*----------------------------------------------------------------------------------------------*/
625: /*@
626:    VecStrideGatherAll - Gathers all the single components from a multi-component vector into
627:    separate vectors.

629:    Collective on Vec

631:    Input Parameter:
632: +  v - the vector
633: -  addv - one of ADD_VALUES,INSERT_VALUES,MAX_VALUES

635:    Output Parameter:
636: .  s - the location where the subvectors are stored

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

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

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

649:    Not optimized; could be easily

651:    Level: advanced

653:    Concepts: gather^into strided vector

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

669:   VecGetLocalSize(v,&n);
670:   VecGetLocalSize(s[0],&n2);
671:   VecGetArrayRead(v,&x);
672:   VecGetBlockSize(v,&bs);
673:   if (bs < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Input vector does not have a valid blocksize set");

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

688:   n =  n/bs;

690:   jj = 0;
691:   if (addv == INSERT_VALUES) {
692:     for (j=0; j<nv; j++) {
693:       for (k=0; k<bss[j]; k++) {
694:         for (i=0; i<n; i++) y[j][i*bss[j] + k] = x[bs*i+jj+k];
695:       }
696:       jj += bss[j];
697:     }
698:   } else if (addv == ADD_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] += x[bs*i+jj+k];
702:       }
703:       jj += bss[j];
704:     }
705: #if !defined(PETSC_USE_COMPLEX)
706:   } else if (addv == MAX_VALUES) {
707:     for (j=0; j<nv; j++) {
708:       for (k=0; k<bss[j]; k++) {
709:         for (i=0; i<n; i++) y[j][i*bss[j] + k] = PetscMax(y[j][i*bss[j] + k],x[bs*i+jj+k]);
710:       }
711:       jj += bss[j];
712:     }
713: #endif
714:   } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown insert type");

716:   VecRestoreArrayRead(v,&x);
717:   for (i=0; i<nv; i++) {
718:     VecRestoreArray(s[i],&y[i]);
719:   }

721:   PetscFree2(y,bss);
722:   return(0);
723: }

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

731:    Collective on Vec

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

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

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

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

747:    Not optimized; could be easily

749:    Level: advanced

751:    Concepts:  scatter^into strided vector

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

766:   VecGetLocalSize(v,&n);
767:   VecGetLocalSize(s[0],&n2);
768:   VecGetArray(v,&x);
769:   VecGetBlockSize(v,&bs);
770:   if (bs < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Input vector does not have a valid blocksize set");

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

785:   n =  n/bs;

787:   jj = 0;
788:   if (addv == INSERT_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] = y[j][i*bss[j] + k];
792:       }
793:       jj += bss[j];
794:     }
795:   } else if (addv == ADD_VALUES) {
796:     for (j=0; j<nv; j++) {
797:       for (k=0; k<bss[j]; k++) {
798:         for (i=0; i<n; i++) x[bs*i+jj+k] += y[j][i*bss[j] + k];
799:       }
800:       jj += bss[j];
801:     }
802: #if !defined(PETSC_USE_COMPLEX)
803:   } else if (addv == MAX_VALUES) {
804:     for (j=0; j<nv; j++) {
805:       for (k=0; k<bss[j]; k++) {
806:         for (i=0; i<n; i++) x[bs*i+jj+k] = PetscMax(x[bs*i+jj+k],y[j][i*bss[j] + k]);
807:       }
808:       jj += bss[j];
809:     }
810: #endif
811:   } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown insert type");

813:   VecRestoreArray(v,&x);
814:   for (i=0; i<nv; i++) {
815:     VecRestoreArray(s[i],&y[i]);
816:   }
817:   PetscFree2(y,bss);
818:   return(0);
819: }

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

827:    Collective on Vec

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

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

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

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

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

847:    Not optimized; could be easily

849:    Level: advanced

851:    Concepts: gather^into strided vector

853: .seealso: VecStrideNorm(), VecStrideScatter(), VecStrideMin(), VecStrideMax(), VecStrideGatherAll(),
854:           VecStrideScatterAll()
855: @*/
856: PetscErrorCode  VecStrideGather(Vec v,PetscInt start,Vec s,InsertMode addv)
857: {

863:   if (start < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative start %D",start);
864:   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);
865:   if (!v->ops->stridegather) SETERRQ(PetscObjectComm((PetscObject)s),PETSC_ERR_SUP,"Not implemented for this Vec class");
866:   (*v->ops->stridegather)(v,start,s,addv);
867:   return(0);
868: }

872: /*@
873:    VecStrideScatter - Scatters a single component from a vector into a multi-component vector.

875:    Collective on Vec

877:    Input Parameter:
878: +  s - the single-component vector
879: .  start - starting point of the subvector (defined by a stride)
880: -  addv - one of ADD_VALUES,INSERT_VALUES,MAX_VALUES

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

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

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

892:    Not optimized; could be easily

894:    Level: advanced

896:    Concepts: scatter^into strided vector

898: .seealso: VecStrideNorm(), VecStrideGather(), VecStrideMin(), VecStrideMax(), VecStrideGatherAll(),
899:           VecStrideScatterAll(), VecStrideSubSetScatter(), VecStrideSubSetGather()
900: @*/
901: PetscErrorCode  VecStrideScatter(Vec s,PetscInt start,Vec v,InsertMode addv)
902: {

908:   if (start < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative start %D",start);
909:   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);
910:   if (!v->ops->stridescatter) SETERRQ(PetscObjectComm((PetscObject)s),PETSC_ERR_SUP,"Not implemented for this Vec class");
911:   (*v->ops->stridescatter)(s,start,v,addv);
912:   return(0);
913: }

917: /*@
918:    VecStrideSubSetGather - Gathers a subset of components from a multi-component vector into
919:    another vector.

921:    Collective on Vec

923:    Input Parameter:
924: +  v - the vector
925: .  nidx - the number of indices
926: .  idxv - the indices of the components 0 <= idxv[0] ...idxv[nidx-1] < bs(v), they need not be sorted
927: .  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
928: -  addv - one of ADD_VALUES,INSERT_VALUES,MAX_VALUES

930:    Output Parameter:
931: .  s - the location where the subvector is stored

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


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

940:    Not optimized; could be easily

942:    Level: advanced

944:    Concepts: gather^into strided vector

946: .seealso: VecStrideNorm(), VecStrideScatter(), VecStrideGather(), VecStrideSubSetScatter(), VecStrideMin(), VecStrideMax(), VecStrideGatherAll(),
947:           VecStrideScatterAll()
948: @*/
949: PetscErrorCode  VecStrideSubSetGather(Vec v,PetscInt nidx,const PetscInt idxv[],const PetscInt idxs[],Vec s,InsertMode addv)
950: {

956:   if (nidx == PETSC_DETERMINE) nidx = s->map->bs;
957:   if (!v->ops->stridesubsetgather) SETERRQ(PetscObjectComm((PetscObject)s),PETSC_ERR_SUP,"Not implemented for this Vec class");
958:   (*v->ops->stridesubsetgather)(v,nidx,idxv,idxs,s,addv);
959:   return(0);
960: }

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

967:    Collective on Vec

969:    Input Parameter:
970: +  s - the smaller-component vector
971: .  nidx - the number of indices in idx
972: .  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
973: .  idxv - the indices of the components in the larger-component vector, 0 <= idx[0] ...idx[nidx-1] < bs(v) they need not be sorted
974: -  addv - one of ADD_VALUES,INSERT_VALUES,MAX_VALUES

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

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

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

985:    Not optimized; could be easily

987:    Level: advanced

989:    Concepts: scatter^into strided vector

991: .seealso: VecStrideNorm(), VecStrideGather(), VecStrideGather(), VecStrideSubSetGather(), VecStrideMin(), VecStrideMax(), VecStrideGatherAll(),
992:           VecStrideScatterAll()
993: @*/
994: PetscErrorCode  VecStrideSubSetScatter(Vec s,PetscInt nidx,const PetscInt idxs[],const PetscInt idxv[],Vec v,InsertMode addv)
995: {

1001:   if (nidx == PETSC_DETERMINE) nidx = s->map->bs;
1002:   if (!v->ops->stridesubsetscatter) SETERRQ(PetscObjectComm((PetscObject)s),PETSC_ERR_SUP,"Not implemented for this Vec class");
1003:   (*v->ops->stridesubsetscatter)(s,nidx,idxs,idxv,v,addv);
1004:   return(0);
1005: }

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

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

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

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

1037:   VecRestoreArrayRead(v,&x);
1038:   VecRestoreArray(s,&y);
1039:   return(0);
1040: }

1044: PetscErrorCode  VecStrideScatter_Default(Vec s,PetscInt start,Vec v,InsertMode addv)
1045: {
1046:   PetscErrorCode    ierr;
1047:   PetscInt          i,n,bs,ns;
1048:   PetscScalar       *x;
1049:   const PetscScalar *y;

1052:   VecGetLocalSize(v,&n);
1053:   VecGetLocalSize(s,&ns);
1054:   VecGetArray(v,&x);
1055:   VecGetArrayRead(s,&y);

1057:   bs = v->map->bs;
1058:   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);
1059:   x += start;
1060:   n  =  n/bs;

1062:   if (addv == INSERT_VALUES) {
1063:     for (i=0; i<n; i++) x[bs*i] = y[i];
1064:   } else if (addv == ADD_VALUES) {
1065:     for (i=0; i<n; i++) x[bs*i] += y[i];
1066: #if !defined(PETSC_USE_COMPLEX)
1067:   } else if (addv == MAX_VALUES) {
1068:     for (i=0; i<n; i++) x[bs*i] = PetscMax(y[i],x[bs*i]);
1069: #endif
1070:   } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown insert type");

1072:   VecRestoreArray(v,&x);
1073:   VecRestoreArrayRead(s,&y);
1074:   return(0);
1075: }

1079: PetscErrorCode  VecStrideSubSetGather_Default(Vec v,PetscInt nidx,const PetscInt idxv[],const PetscInt idxs[],Vec s,InsertMode addv)
1080: {
1081:   PetscErrorCode    ierr;
1082:   PetscInt          i,j,n,bs,bss,ns;
1083:   const PetscScalar *x;
1084:   PetscScalar       *y;

1087:   VecGetLocalSize(v,&n);
1088:   VecGetLocalSize(s,&ns);
1089:   VecGetArrayRead(v,&x);
1090:   VecGetArray(s,&y);

1092:   bs  = v->map->bs;
1093:   bss = s->map->bs;
1094:   n  =  n/bs;

1096: #if defined(PETSC_DEBUG)
1097:   if (n != ns/bss) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Incompatible layout of vectors");
1098:   for (j=0; j<nidx; j++) {
1099:     if (idxv[j] < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"idx[%D] %D is negative",j,idxv[j]);
1100:     if (idxv[j] >= bs) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"idx[%D] %D is greater than or equal to vector blocksize %D",j,idxv[j],bs);
1101:   }
1102:   if (!idxs && bss != nidx) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Must provide idxs when not gathering into all locations");
1103: #endif

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

1139:   VecRestoreArrayRead(v,&x);
1140:   VecRestoreArray(s,&y);
1141:   return(0);
1142: }

1146: PetscErrorCode  VecStrideSubSetScatter_Default(Vec s,PetscInt nidx,const PetscInt idxs[],const PetscInt idxv[],Vec v,InsertMode addv)
1147: {
1148:   PetscErrorCode    ierr;
1149:   PetscInt          j,i,n,bs,ns,bss;
1150:   PetscScalar       *x;
1151:   const PetscScalar *y;

1154:   VecGetLocalSize(v,&n);
1155:   VecGetLocalSize(s,&ns);
1156:   VecGetArray(v,&x);
1157:   VecGetArrayRead(s,&y);

1159:   bs  = v->map->bs;
1160:   bss = s->map->bs;
1161:   n  =  n/bs;

1163: #if defined(PETSC_DEBUG)
1164:   if (n != ns/bss) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Incompatible layout of vectors");
1165:   for (j=0; j<bss; j++) {
1166:     if (idx[j] < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"idx[%D] %D is negative",j,idx[j]);
1167:     if (idx[j] >= bs) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"idx[%D] %D is greater than or equal to vector blocksize %D",j,idx[j],bs);
1168:   }
1169:   if (!idxs && bss != nidx) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Must provide idxs when not scattering from all locations");
1170: #endif

1172:   if (addv == INSERT_VALUES) {
1173:     if (!idxs) {
1174:       for (i=0; i<n; i++) {
1175:         for (j=0; j<bss; j++) x[bs*i + idxv[j]] = y[bss*i+j];
1176:       }
1177:     } else {
1178:       for (i=0; i<n; i++) {
1179:         for (j=0; j<bss; j++) x[bs*i + idxv[j]] = y[bss*i+idxs[j]];
1180:       }
1181:     }
1182:   } else if (addv == ADD_VALUES) {
1183:     if (!idxs) {
1184:       for (i=0; i<n; i++) {
1185:         for (j=0; j<bss; j++) x[bs*i + idxv[j]] += y[bss*i+j];
1186:       }
1187:     } else {
1188:       for (i=0; i<n; i++) {
1189:         for (j=0; j<bss; j++) x[bs*i + idxv[j]] += y[bss*i+idxs[j]];
1190:       }
1191:     }
1192: #if !defined(PETSC_USE_COMPLEX)
1193:   } else if (addv == MAX_VALUES) {
1194:     if (!idxs) {
1195:       for (i=0; i<n; i++) {
1196:         for (j=0; j<bss; j++) x[bs*i + idxv[j]] = PetscMax(y[bss*i+j],x[bs*i + idxv[j]]);
1197:       }
1198:     } else {
1199:       for (i=0; i<n; i++) {
1200:         for (j=0; j<bss; j++) x[bs*i + idxv[j]] = PetscMax(y[bss*i+idxs[j]],x[bs*i + idxv[j]]);
1201:       }
1202:     }
1203: #endif
1204:   } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown insert type");

1206:   VecRestoreArray(v,&x);
1207:   VecRestoreArrayRead(s,&y);
1208:   return(0);
1209: }

1213: PetscErrorCode VecReciprocal_Default(Vec v)
1214: {
1216:   PetscInt       i,n;
1217:   PetscScalar    *x;

1220:   VecGetLocalSize(v,&n);
1221:   VecGetArray(v,&x);
1222:   for (i=0; i<n; i++) {
1223:     if (x[i] != (PetscScalar)0.0) x[i] = (PetscScalar)1.0/x[i];
1224:   }
1225:   VecRestoreArray(v,&x);
1226:   return(0);
1227: }

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

1234:   Not collective

1236:   Input Parameter:
1237: . v - The vector

1239:   Output Parameter:
1240: . v - The vector of exponents

1242:   Level: beginner

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

1246: .keywords: vector, sqrt, square root
1247: @*/
1248: PetscErrorCode  VecExp(Vec v)
1249: {
1250:   PetscScalar    *x;
1251:   PetscInt       i, n;

1256:   if (v->ops->exp) {
1257:     (*v->ops->exp)(v);
1258:   } else {
1259:     VecGetLocalSize(v, &n);
1260:     VecGetArray(v, &x);
1261:     for (i = 0; i < n; i++) x[i] = PetscExpScalar(x[i]);
1262:     VecRestoreArray(v, &x);
1263:   }
1264:   return(0);
1265: }

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

1272:   Not collective

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

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

1280:   Level: beginner

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

1284: .keywords: vector, sqrt, square root
1285: @*/
1286: PetscErrorCode  VecLog(Vec v)
1287: {
1288:   PetscScalar    *x;
1289:   PetscInt       i, n;

1294:   if (v->ops->log) {
1295:     (*v->ops->log)(v);
1296:   } else {
1297:     VecGetLocalSize(v, &n);
1298:     VecGetArray(v, &x);
1299:     for (i = 0; i < n; i++) x[i] = PetscLogScalar(x[i]);
1300:     VecRestoreArray(v, &x);
1301:   }
1302:   return(0);
1303: }

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

1310:   Not collective

1312:   Input Parameter:
1313: . v - The vector

1315:   Output Parameter:
1316: . v - The vector square root

1318:   Level: beginner

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

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

1324: .keywords: vector, sqrt, square root
1325: @*/
1326: PetscErrorCode  VecSqrtAbs(Vec v)
1327: {
1328:   PetscScalar    *x;
1329:   PetscInt       i, n;

1334:   if (v->ops->sqrt) {
1335:     (*v->ops->sqrt)(v);
1336:   } else {
1337:     VecGetLocalSize(v, &n);
1338:     VecGetArray(v, &x);
1339:     for (i = 0; i < n; i++) x[i] = PetscSqrtReal(PetscAbsScalar(x[i]));
1340:     VecRestoreArray(v, &x);
1341:   }
1342:   return(0);
1343: }

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

1350:   Collective on Vec

1352:   Input Parameter:
1353: + s - first vector
1354: - t - second vector

1356:   Output Parameter:
1357: + dp - s'conj(t)
1358: - nm - t'conj(t)

1360:   Level: advanced

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


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

1367: .keywords: vector, sqrt, square root
1368: @*/
1369: PetscErrorCode  VecDotNorm2(Vec s,Vec t,PetscScalar *dp, PetscReal *nm)
1370: {
1371:   const PetscScalar *sx, *tx;
1372:   PetscScalar       dpx = 0.0, nmx = 0.0,work[2],sum[2];
1373:   PetscInt          i, n;
1374:   PetscErrorCode    ierr;

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

1387:   PetscLogEventBarrierBegin(VEC_DotNormBarrier,s,t,0,0,PetscObjectComm((PetscObject)s));
1388:   if (s->ops->dotnorm2) {
1389:     (*s->ops->dotnorm2)(s,t,dp,&dpx);
1390:     *nm  = PetscRealPart(dpx);
1391:   } else {
1392:     VecGetLocalSize(s, &n);
1393:     VecGetArrayRead(s, &sx);
1394:     VecGetArrayRead(t, &tx);

1396:     for (i = 0; i<n; i++) {
1397:       dpx += sx[i]*PetscConj(tx[i]);
1398:       nmx += tx[i]*PetscConj(tx[i]);
1399:     }
1400:     work[0] = dpx;
1401:     work[1] = nmx;

1403:     MPI_Allreduce(work,sum,2,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)s));
1404:     *dp  = sum[0];
1405:     *nm  = PetscRealPart(sum[1]);

1407:     VecRestoreArrayRead(t, &tx);
1408:     VecRestoreArrayRead(s, &sx);
1409:     PetscLogFlops(4.0*n);
1410:   }
1411:   PetscLogEventBarrierEnd(VEC_DotNormBarrier,s,t,0,0,PetscObjectComm((PetscObject)s));
1412:   return(0);
1413: }

1417: /*@
1418:    VecSum - Computes the sum of all the components of a vector.

1420:    Collective on Vec

1422:    Input Parameter:
1423: .  v - the vector

1425:    Output Parameter:
1426: .  sum - the result

1428:    Level: beginner

1430:    Concepts: sum^of vector entries

1432: .seealso: VecNorm()
1433: @*/
1434: PetscErrorCode  VecSum(Vec v,PetscScalar *sum)
1435: {
1436:   PetscErrorCode    ierr;
1437:   PetscInt          i,n;
1438:   const PetscScalar *x;
1439:   PetscScalar       lsum = 0.0;

1444:   VecGetLocalSize(v,&n);
1445:   VecGetArrayRead(v,&x);
1446:   for (i=0; i<n; i++) lsum += x[i];
1447:   MPI_Allreduce(&lsum,sum,1,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)v));
1448:   VecRestoreArrayRead(v,&x);
1449:   return(0);
1450: }

1454: /*@
1455:    VecShift - Shifts all of the components of a vector by computing
1456:    x[i] = x[i] + shift.

1458:    Logically Collective on Vec

1460:    Input Parameters:
1461: +  v - the vector
1462: -  shift - the shift

1464:    Output Parameter:
1465: .  v - the shifted vector

1467:    Level: intermediate

1469:    Concepts: vector^adding constant

1471: @*/
1472: PetscErrorCode  VecShift(Vec v,PetscScalar shift)
1473: {
1475:   PetscInt       i,n;
1476:   PetscScalar    *x;

1481:   VecLocked(v,1);

1483:   if (v->ops->shift) {
1484:     (*v->ops->shift)(v);
1485:   } else {
1486:     VecGetLocalSize(v,&n);
1487:     VecGetArray(v,&x);
1488:     for (i=0; i<n; i++) x[i] += shift;
1489:     VecRestoreArray(v,&x);
1490:   }
1491:   return(0);
1492: }

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

1499:    Logically Collective on Vec

1501:    Input Parameters:
1502: .  v - the vector

1504:    Level: intermediate

1506:    Concepts: vector^absolute value

1508: @*/
1509: PetscErrorCode  VecAbs(Vec v)
1510: {
1512:   PetscInt       i,n;
1513:   PetscScalar    *x;

1517:   VecLocked(v,1);
1518: 
1519:   if (v->ops->abs) {
1520:     (*v->ops->abs)(v);
1521:   } else {
1522:     VecGetLocalSize(v,&n);
1523:     VecGetArray(v,&x);
1524:     for (i=0; i<n; i++) x[i] = PetscAbsScalar(x[i]);
1525:     VecRestoreArray(v,&x);
1526:   }
1527:   return(0);
1528: }

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

1535:   Input Parameters:
1536: + vec   - The vector
1537: . order - The ordering
1538: - inv   - The flag for inverting the permutation

1540:   Level: beginner

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

1544: .seealso: MatPermute()
1545: .keywords: vec, permute
1546: @*/
1547: PetscErrorCode  VecPermute(Vec x, IS row, PetscBool inv)
1548: {
1549:   PetscScalar    *array, *newArray;
1550:   const PetscInt *idx;
1551:   PetscInt       i,rstart,rend;

1555:   VecLocked(x,1);
1556: 
1557:   VecGetOwnershipRange(x,&rstart,&rend);
1558:   ISGetIndices(row, &idx);
1559:   VecGetArray(x, &array);
1560:   PetscMalloc1(x->map->n, &newArray);
1561: #if defined(PETSC_USE_DEBUG)
1562:   for (i = 0; i < x->map->n; i++) {
1563:     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]);
1564:   }
1565: #endif
1566:   if (!inv) {
1567:     for (i = 0; i < x->map->n; i++) newArray[i] = array[idx[i]-rstart];
1568:   } else {
1569:     for (i = 0; i < x->map->n; i++) newArray[idx[i]-rstart] = array[i];
1570:   }
1571:   VecRestoreArray(x, &array);
1572:   ISRestoreIndices(row, &idx);
1573:   VecReplaceArray(x, newArray);
1574:   return(0);
1575: }

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

1584:    Collective on Vec

1586:    Input Parameters:
1587: +  vec1 - the first vector
1588: -  vec2 - the second vector

1590:    Output Parameter:
1591: .  flg - PETSC_TRUE if the vectors are equal; PETSC_FALSE otherwise.

1593:    Level: intermediate

1595:    Concepts: equal^two vectors
1596:    Concepts: vector^equality

1598: @*/
1599: PetscErrorCode  VecEqual(Vec vec1,Vec vec2,PetscBool  *flg)
1600: {
1601:   const PetscScalar  *v1,*v2;
1602:   PetscErrorCode     ierr;
1603:   PetscInt           n1,n2,N1,N2;
1604:   PetscBool          flg1;

1610:   if (vec1 == vec2) *flg = PETSC_TRUE;
1611:   else {
1612:     VecGetSize(vec1,&N1);
1613:     VecGetSize(vec2,&N2);
1614:     if (N1 != N2) flg1 = PETSC_FALSE;
1615:     else {
1616:       VecGetLocalSize(vec1,&n1);
1617:       VecGetLocalSize(vec2,&n2);
1618:       if (n1 != n2) flg1 = PETSC_FALSE;
1619:       else {
1620:         VecGetArrayRead(vec1,&v1);
1621:         VecGetArrayRead(vec2,&v2);
1622:         PetscMemcmp(v1,v2,n1*sizeof(PetscScalar),&flg1);
1623:         VecRestoreArrayRead(vec1,&v1);
1624:         VecRestoreArrayRead(vec2,&v2);
1625:       }
1626:     }
1627:     /* combine results from all processors */
1628:     MPI_Allreduce(&flg1,flg,1,MPIU_BOOL,MPI_MIN,PetscObjectComm((PetscObject)vec1));
1629:   }
1630:   return(0);
1631: }

1635: /*@
1636:    VecUniqueEntries - Compute the number of unique entries, and those entries

1638:    Collective on Vec

1640:    Input Parameter:
1641: .  vec - the vector

1643:    Output Parameters:
1644: +  n - The number of unique entries
1645: -  e - The entries

1647:    Level: intermediate

1649: @*/
1650: PetscErrorCode  VecUniqueEntries(Vec vec, PetscInt *n, PetscScalar **e)
1651: {
1652:   const PetscScalar *v;
1653:   PetscScalar       *tmp, *vals;
1654:   PetscMPIInt       *N, *displs, l;
1655:   PetscInt          ng, m, i, j, p;
1656:   PetscMPIInt       size;
1657:   PetscErrorCode    ierr;

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