Actual source code: vinv.c

petsc-3.5.1 2014-08-06
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 multiply each subvector entry by

 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;


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

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

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

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

 66:    Logically Collective on Vec

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

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

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

 79:    Level: advanced

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

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


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

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

107:   VecRestoreArray(v,&x);
108:   return(0);
109: }

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

117:    Collective on Vec

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

124:    Output Parameter:
125: .  norm - the norm

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

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

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

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

139:    Level: advanced

141:    Concepts: norm^on stride of vector
142:    Concepts: stride^norm

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

157:   VecGetLocalSize(v,&n);
158:   VecGetArray(v,&x);
159:   PetscObjectGetComm((PetscObject)v,&comm);

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

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

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

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

197:    Collective on Vec

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

203:    Output Parameter:
204: +  index - the location where the maximum occurred  (pass NULL if not required)
205: -  nrm - the max

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

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

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

218:    Level: advanced

220:    Concepts: maximum^on stride of vector
221:    Concepts: stride^maximum

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


237:   VecGetLocalSize(v,&n);
238:   VecGetArray(v,&x);
239:   PetscObjectGetComm((PetscObject)v,&comm);

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

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

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

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

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

279:    Collective on Vec

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

285:    Output Parameter:
286: +  idex - the location where the minimum occurred. (pass NULL if not required)
287: -  nrm - the min

289:    Level: advanced

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

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

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

302:    Concepts: minimum^on stride of vector
303:    Concepts: stride^minimum

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


319:   VecGetLocalSize(v,&n);
320:   VecGetArray(v,&x);
321:   PetscObjectGetComm((PetscObject)v,&comm);

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

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

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

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

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

361:    Logically Collective on Vec

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

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


372:    Level: advanced

374:    Concepts: scale^on stride of vector
375:    Concepts: stride^scale

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

388:   VecGetLocalSize(v,&n);
389:   VecGetArray(v,&x);
390:   VecGetBlockSize(v,&bs);

392:   /* need to provide optimized code for each bs */
393:   for (i=0; i<n; i+=bs) {
394:     for (j=0; j<bs; j++) x[i+j] *= scales[j];
395:   }
396:   VecRestoreArray(v,&x);
397:   return(0);
398: }

402: /*@
403:    VecStrideNormAll - Computes the norms  subvectors of a vector defined
404:    by a starting point and a stride.

406:    Collective on Vec

408:    Input Parameter:
409: +  v - the vector
410: -  ntype - type of norm, one of NORM_1, NORM_2, NORM_INFINITY

412:    Output Parameter:
413: .  nrm - the norms

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

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

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

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

427:    Level: advanced

429:    Concepts: norm^on stride of vector
430:    Concepts: stride^norm

432: .seealso: VecNorm(), VecStrideGather(), VecStrideScatter(), VecStrideMin(), VecStrideMax()
433: @*/
434: PetscErrorCode  VecStrideNormAll(Vec v,NormType ntype,PetscReal nrm[])
435: {
437:   PetscInt       i,j,n,bs;
438:   PetscScalar    *x;
439:   PetscReal      tnorm[128];
440:   MPI_Comm       comm;

445:   VecGetLocalSize(v,&n);
446:   VecGetArray(v,&x);
447:   PetscObjectGetComm((PetscObject)v,&comm);

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

452:   if (ntype == NORM_2) {
453:     PetscScalar sum[128];
454:     for (j=0; j<bs; j++) sum[j] = 0.0;
455:     for (i=0; i<n; i+=bs) {
456:       for (j=0; j<bs; j++) sum[j] += x[i+j]*(PetscConj(x[i+j]));
457:     }
458:     for (j=0; j<bs; j++) tnorm[j]  = PetscRealPart(sum[j]);

460:     MPI_Allreduce(tnorm,nrm,bs,MPIU_REAL,MPIU_SUM,comm);
461:     for (j=0; j<bs; j++) nrm[j] = PetscSqrtReal(nrm[j]);
462:   } else if (ntype == NORM_1) {
463:     for (j=0; j<bs; j++) tnorm[j] = 0.0;

465:     for (i=0; i<n; i+=bs) {
466:       for (j=0; j<bs; j++) tnorm[j] += PetscAbsScalar(x[i+j]);
467:     }

469:     MPI_Allreduce(tnorm,nrm,bs,MPIU_REAL,MPIU_SUM,comm);
470:   } else if (ntype == NORM_INFINITY) {
471:     PetscReal tmp;
472:     for (j=0; j<bs; j++) tnorm[j] = 0.0;

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

489: /*@
490:    VecStrideMaxAll - Computes the maximums of subvectors of a vector defined
491:    by a starting point and a stride and optionally its location.

493:    Collective on Vec

495:    Input Parameter:
496: .  v - the vector

498:    Output Parameter:
499: +  index - the location where the maximum occurred (not supported, pass NULL,
500:            if you need this, send mail to petsc-maint@mcs.anl.gov to request it)
501: -  nrm - the maximums

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

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

511:    Level: advanced

513:    Concepts: maximum^on stride of vector
514:    Concepts: stride^maximum

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

529:   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");
530:   VecGetLocalSize(v,&n);
531:   VecGetArray(v,&x);
532:   PetscObjectGetComm((PetscObject)v,&comm);

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

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

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

550:   VecRestoreArray(v,&x);
551:   return(0);
552: }

556: /*@
557:    VecStrideMinAll - Computes the minimum of subvector of a vector defined
558:    by a starting point and a stride and optionally its location.

560:    Collective on Vec

562:    Input Parameter:
563: .  v - the vector

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

570:    Level: advanced

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

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

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

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

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

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

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

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

617:   VecRestoreArray(v,&x);
618:   return(0);
619: }

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

628:    Collective on Vec

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

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

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

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

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

648:    Not optimized; could be easily

650:    Level: advanced

652:    Concepts: gather^into strided vector

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

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

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

686:   n =  n/bs;

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

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

719:   PetscFree2(y,bss);
720:   return(0);
721: }

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

729:    Collective on Vec

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

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

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

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

745:    Not optimized; could be easily

747:    Level: advanced

749:    Concepts:  scatter^into strided vector

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

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

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

783:   n =  n/bs;

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

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

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

825:    Collective on Vec

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

832:    Output Parameter:
833: .  s - the location where the subvector is stored

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

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

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

845:    Not optimized; could be easily

847:    Level: advanced

849:    Concepts: gather^into strided vector

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

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

870: /*@
871:    VecStrideScatter - Scatters a single component from a vector into a multi-component vector.

873:    Collective on Vec

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

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

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

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

890:    Not optimized; could be easily

892:    Level: advanced

894:    Concepts: scatter^into strided vector

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

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

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

919:    Collective on Vec

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

928:    Output Parameter:
929: .  s - the location where the subvector is stored

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


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

938:    Not optimized; could be easily

940:    Level: advanced

942:    Concepts: gather^into strided vector

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

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

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

965:    Collective on Vec

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

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

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

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

983:    Not optimized; could be easily

985:    Level: advanced

987:    Concepts: scatter^into strided vector

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

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

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

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

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

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

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

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

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

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

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

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

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

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

1090:   bs  = v->map->bs;
1091:   bss = s->map->bs;
1092:   n  =  n/bs;

1094: #if defined(PETSC_DEBUG)
1095:   if (n != ns/bss) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Incompatible layout of vectors");
1096:   for (j=0; j<nidx; j++) {
1097:     if (idxv[j] < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"idx[%D] %D is negative",j,idxv[j]);
1098:     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);
1099:   }
1100:   if (!idxs && bss != nidx) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Must provide idxs when not gathering into all locations");
1101: #endif

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

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

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

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

1157:   bs  = v->map->bs;
1158:   bss = s->map->bs;
1159:   n  =  n/bs;

1161: #if defined(PETSC_DEBUG)
1162:   if (n != ns/bss) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Incompatible layout of vectors");
1163:   for (j=0; j<bss; j++) {
1164:     if (idx[j] < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"idx[%D] %D is negative",j,idx[j]);
1165:     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);
1166:   }
1167:   if (!idxs && bss != nidx) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Must provide idxs when not scattering from all locations");
1168: #endif

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

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

1211: PetscErrorCode VecReciprocal_Default(Vec v)
1212: {
1214:   PetscInt       i,n;
1215:   PetscScalar    *x;

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

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

1232:   Not collective

1234:   Input Parameter:
1235: . v - The vector

1237:   Output Parameter:
1238: . v - The vector of exponents

1240:   Level: beginner

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

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

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

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

1270:   Not collective

1272:   Input Parameter:
1273: . v - The vector

1275:   Output Parameter:
1276: . v - The vector of logs

1278:   Level: beginner

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

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

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

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

1308:   Not collective

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

1313:   Output Parameter:
1314: . v - The vector square root

1316:   Level: beginner

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

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

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

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

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

1348:   Collective on Vec

1350:   Input Parameter:
1351: + s - first vector
1352: - t - second vector

1354:   Output Parameter:
1355: + dp - s'conj(t)
1356: - nm - t'conj(t)

1358:   Level: advanced

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


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

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

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

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

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

1400:     MPI_Allreduce(work,sum,2,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)s));
1401:     *dp  = sum[0];
1402:     *nm  = PetscRealPart(sum[1]);

1404:     VecRestoreArray(t, &tx);
1405:     VecRestoreArray(s, &sx);
1406:     PetscLogFlops(4.0*n);
1407:   }
1408:   PetscLogEventBarrierEnd(VEC_DotNormBarrier,s,t,0,0,PetscObjectComm((PetscObject)s));
1409:   return(0);
1410: }

1414: /*@
1415:    VecSum - Computes the sum of all the components of a vector.

1417:    Collective on Vec

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

1422:    Output Parameter:
1423: .  sum - the result

1425:    Level: beginner

1427:    Concepts: sum^of vector entries

1429: .seealso: VecNorm()
1430: @*/
1431: PetscErrorCode  VecSum(Vec v,PetscScalar *sum)
1432: {
1434:   PetscInt       i,n;
1435:   PetscScalar    *x,lsum = 0.0;

1440:   VecGetLocalSize(v,&n);
1441:   VecGetArray(v,&x);
1442:   for (i=0; i<n; i++) lsum += x[i];
1443:   MPI_Allreduce(&lsum,sum,1,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)v));
1444:   VecRestoreArray(v,&x);
1445:   return(0);
1446: }

1450: /*@
1451:    VecShift - Shifts all of the components of a vector by computing
1452:    x[i] = x[i] + shift.

1454:    Logically Collective on Vec

1456:    Input Parameters:
1457: +  v - the vector
1458: -  shift - the shift

1460:    Output Parameter:
1461: .  v - the shifted vector

1463:    Level: intermediate

1465:    Concepts: vector^adding constant

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

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

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

1493:    Logically Collective on Vec

1495:    Input Parameters:
1496: .  v - the vector

1498:    Level: intermediate

1500:    Concepts: vector^absolute value

1502: @*/
1503: PetscErrorCode  VecAbs(Vec v)
1504: {
1506:   PetscInt       i,n;
1507:   PetscScalar    *x;

1511:   if (v->ops->abs) {
1512:     (*v->ops->abs)(v);
1513:   } else {
1514:     VecGetLocalSize(v,&n);
1515:     VecGetArray(v,&x);
1516:     for (i=0; i<n; i++) x[i] = PetscAbsScalar(x[i]);
1517:     VecRestoreArray(v,&x);
1518:   }
1519:   return(0);
1520: }

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

1527:   Input Parameters:
1528: + vec   - The vector
1529: . order - The ordering
1530: - inv   - The flag for inverting the permutation

1532:   Level: beginner

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

1536: .seealso: MatPermute()
1537: .keywords: vec, permute
1538: @*/
1539: PetscErrorCode  VecPermute(Vec x, IS row, PetscBool inv)
1540: {
1541:   PetscScalar    *array, *newArray;
1542:   const PetscInt *idx;
1543:   PetscInt       i,rstart,rend;

1547:   VecGetOwnershipRange(x,&rstart,&rend);
1548:   ISGetIndices(row, &idx);
1549:   VecGetArray(x, &array);
1550:   PetscMalloc1(x->map->n, &newArray);
1551: #if defined(PETSC_USE_DEBUG)
1552:   for (i = 0; i < x->map->n; i++) {
1553:     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]);
1554:   }
1555: #endif
1556:   if (!inv) {
1557:     for (i = 0; i < x->map->n; i++) newArray[i] = array[idx[i]-rstart];
1558:   } else {
1559:     for (i = 0; i < x->map->n; i++) newArray[idx[i]-rstart] = array[i];
1560:   }
1561:   VecRestoreArray(x, &array);
1562:   ISRestoreIndices(row, &idx);
1563:   VecReplaceArray(x, newArray);
1564:   return(0);
1565: }

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

1574:    Collective on Vec

1576:    Input Parameters:
1577: +  vec1 - the first vector
1578: -  vec2 - the second vector

1580:    Output Parameter:
1581: .  flg - PETSC_TRUE if the vectors are equal; PETSC_FALSE otherwise.

1583:    Level: intermediate

1585:    Concepts: equal^two vectors
1586:    Concepts: vector^equality

1588: @*/
1589: PetscErrorCode  VecEqual(Vec vec1,Vec vec2,PetscBool  *flg)
1590: {
1591:   const PetscScalar  *v1,*v2;
1592:   PetscErrorCode     ierr;
1593:   PetscInt           n1,n2,N1,N2;
1594:   PetscBool          flg1;

1600:   if (vec1 == vec2) *flg = PETSC_TRUE;
1601:   else {
1602:     VecGetSize(vec1,&N1);
1603:     VecGetSize(vec2,&N2);
1604:     if (N1 != N2) flg1 = PETSC_FALSE;
1605:     else {
1606:       VecGetLocalSize(vec1,&n1);
1607:       VecGetLocalSize(vec2,&n2);
1608:       if (n1 != n2) flg1 = PETSC_FALSE;
1609:       else {
1610:         VecGetArrayRead(vec1,&v1);
1611:         VecGetArrayRead(vec2,&v2);
1612:         PetscMemcmp(v1,v2,n1*sizeof(PetscScalar),&flg1);
1613:         VecRestoreArrayRead(vec1,&v1);
1614:         VecRestoreArrayRead(vec2,&v2);
1615:       }
1616:     }
1617:     /* combine results from all processors */
1618:     MPI_Allreduce(&flg1,flg,1,MPIU_BOOL,MPI_MIN,PetscObjectComm((PetscObject)vec1));
1619:   }
1620:   return(0);
1621: }

1625: /*@
1626:    VecUniqueEntries - Compute the number of unique entries, and those entries

1628:    Collective on Vec

1630:    Input Parameter:
1631: .  vec - the vector

1633:    Output Parameters:
1634: +  n - The number of unique entries
1635: -  e - The entries

1637:    Level: intermediate

1639: @*/
1640: PetscErrorCode  VecUniqueEntries(Vec vec, PetscInt *n, PetscScalar **e)
1641: {
1642:   PetscScalar   *v, *tmp, *vals;
1643:   PetscMPIInt   *N, *displs, l;
1644:   PetscInt       ng, m, i, j, p;
1645:   PetscMPIInt    size;

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