Actual source code: vinv.c

petsc-3.8.3 2017-12-09
Report Typos and Errors

  2: /*
  3:      Some useful vector utility functions.
  4: */
  5:  #include <../src/vec/vec/impls/mpi/pvecimpl.h>
  6: extern MPI_Op MPIU_MAXINDEX_OP;
  7: extern MPI_Op MPIU_MININDEX_OP;

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

 13:    Logically Collective on Vec

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

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

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

 26:    Level: advanced

 28:    Concepts: scale^on stride of vector
 29:    Concepts: stride^scale

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

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

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

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

 59: /*@
 60:    VecStrideScale - Scales a subvector of a vector defined
 61:    by a starting point and a stride.

 63:    Logically Collective on Vec

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

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

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

 76:    Level: advanced

 78:    Concepts: scale^on stride of vector
 79:    Concepts: stride^scale

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

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

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

105:   VecRestoreArray(v,&x);
106:   return(0);
107: }

109: /*@
110:    VecStrideNorm - Computes the norm of subvector of a vector defined
111:    by a starting point and a stride.

113:    Collective on Vec

115:    Input Parameter:
116: +  v - the vector
117: .  start - starting point of the subvector (defined by a stride)
118: -  ntype - type of norm, one of NORM_1, NORM_2, NORM_INFINITY

120:    Output Parameter:
121: .  norm - the norm

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

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

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

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

135:    Level: advanced

137:    Concepts: norm^on stride of vector
138:    Concepts: stride^norm

140: .seealso: VecNorm(), VecStrideGather(), VecStrideScatter(), VecStrideMin(), VecStrideMax()
141: @*/
142: PetscErrorCode  VecStrideNorm(Vec v,PetscInt start,NormType ntype,PetscReal *nrm)
143: {
144:   PetscErrorCode    ierr;
145:   PetscInt          i,n,bs;
146:   const PetscScalar *x;
147:   PetscReal         tnorm;
148:   MPI_Comm          comm;

153:   VecGetLocalSize(v,&n);
154:   VecGetArrayRead(v,&x);
155:   PetscObjectGetComm((PetscObject)v,&comm);

157:   VecGetBlockSize(v,&bs);
158:   if (start < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative start %D",start);
159:   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);
160:   x += start;

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

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

187: /*@
188:    VecStrideMax - Computes the maximum of subvector of a vector defined
189:    by a starting point and a stride and optionally its location.

191:    Collective on Vec

193:    Input Parameter:
194: +  v - the vector
195: -  start - starting point of the subvector (defined by a stride)

197:    Output Parameter:
198: +  index - the location where the maximum occurred  (pass NULL if not required)
199: -  nrm - the maximum value in the subvector

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

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

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

212:    Level: advanced

214:    Concepts: maximum^on stride of vector
215:    Concepts: stride^maximum

217: .seealso: VecMax(), VecStrideNorm(), VecStrideGather(), VecStrideScatter(), VecStrideMin()
218: @*/
219: PetscErrorCode  VecStrideMax(Vec v,PetscInt start,PetscInt *idex,PetscReal *nrm)
220: {
221:   PetscErrorCode    ierr;
222:   PetscInt          i,n,bs,id;
223:   const PetscScalar *x;
224:   PetscReal         max,tmp;
225:   MPI_Comm          comm;


231:   VecGetLocalSize(v,&n);
232:   VecGetArrayRead(v,&x);
233:   PetscObjectGetComm((PetscObject)v,&comm);

235:   VecGetBlockSize(v,&bs);
236:   if (start < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative start %D",start);
237:   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);
238:   x += start;

240:   id = -1;
241:   if (!n) max = PETSC_MIN_REAL;
242:   else {
243:     id  = 0;
244:     max = PetscRealPart(x[0]);
245:     for (i=bs; i<n; i+=bs) {
246:       if ((tmp = PetscRealPart(x[i])) > max) { max = tmp; id = i;}
247:     }
248:   }
249:   VecRestoreArrayRead(v,&x);

251:   if (!idex) {
252:     MPIU_Allreduce(&max,nrm,1,MPIU_REAL,MPIU_MAX,comm);
253:   } else {
254:     PetscReal in[2],out[2];
255:     PetscInt  rstart;

257:     VecGetOwnershipRange(v,&rstart,NULL);
258:     in[0] = max;
259:     in[1] = rstart+id+start;
260:     MPIU_Allreduce(in,out,2,MPIU_REAL,MPIU_MAXINDEX_OP,PetscObjectComm((PetscObject)v));
261:     *nrm  = out[0];
262:     *idex = (PetscInt)out[1];
263:   }
264:   return(0);
265: }

267: /*@
268:    VecStrideMin - Computes the minimum of subvector of a vector defined
269:    by a starting point and a stride and optionally its location.

271:    Collective on Vec

273:    Input Parameter:
274: +  v - the vector
275: -  start - starting point of the subvector (defined by a stride)

277:    Output Parameter:
278: +  idex - the location where the minimum occurred. (pass NULL if not required)
279: -  nrm - the minimum value in the subvector

281:    Level: advanced

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

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

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

294:    Concepts: minimum^on stride of vector
295:    Concepts: stride^minimum

297: .seealso: VecMin(), VecStrideNorm(), VecStrideGather(), VecStrideScatter(), VecStrideMax()
298: @*/
299: PetscErrorCode  VecStrideMin(Vec v,PetscInt start,PetscInt *idex,PetscReal *nrm)
300: {
301:   PetscErrorCode    ierr;
302:   PetscInt          i,n,bs,id;
303:   const PetscScalar *x;
304:   PetscReal         min,tmp;
305:   MPI_Comm          comm;


311:   VecGetLocalSize(v,&n);
312:   VecGetArrayRead(v,&x);
313:   PetscObjectGetComm((PetscObject)v,&comm);

315:   VecGetBlockSize(v,&bs);
316:   if (start < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative start %D",start);
317:   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);
318:   x += start;

320:   id = -1;
321:   if (!n) min = PETSC_MAX_REAL;
322:   else {
323:     id = 0;
324:     min = PetscRealPart(x[0]);
325:     for (i=bs; i<n; i+=bs) {
326:       if ((tmp = PetscRealPart(x[i])) < min) { min = tmp; id = i;}
327:     }
328:   }
329:   VecRestoreArrayRead(v,&x);

331:   if (!idex) {
332:     MPIU_Allreduce(&min,nrm,1,MPIU_REAL,MPIU_MIN,comm);
333:   } else {
334:     PetscReal in[2],out[2];
335:     PetscInt  rstart;

337:     VecGetOwnershipRange(v,&rstart,NULL);
338:     in[0] = min;
339:     in[1] = rstart+id;
340:     MPIU_Allreduce(in,out,2,MPIU_REAL,MPIU_MININDEX_OP,PetscObjectComm((PetscObject)v));
341:     *nrm  = out[0];
342:     *idex = (PetscInt)out[1];
343:   }
344:   return(0);
345: }

347: /*@
348:    VecStrideScaleAll - Scales the subvectors of a vector defined
349:    by a starting point and a stride.

351:    Logically Collective on Vec

353:    Input Parameter:
354: +  v - the vector
355: -  scales - values to multiply each subvector entry by

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

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


364:    Level: advanced

366:    Concepts: scale^on stride of vector
367:    Concepts: stride^scale

369: .seealso: VecNorm(), VecStrideScale(), VecScale(), VecStrideGather(), VecStrideScatter(), VecStrideMin(), VecStrideMax()
370: @*/
371: PetscErrorCode  VecStrideScaleAll(Vec v,const PetscScalar *scales)
372: {
374:   PetscInt       i,j,n,bs;
375:   PetscScalar    *x;

380:   VecLocked(v,1);
381: 
382:   VecGetLocalSize(v,&n);
383:   VecGetArray(v,&x);
384:   VecGetBlockSize(v,&bs);

386:   /* need to provide optimized code for each bs */
387:   for (i=0; i<n; i+=bs) {
388:     for (j=0; j<bs; j++) x[i+j] *= scales[j];
389:   }
390:   VecRestoreArray(v,&x);
391:   return(0);
392: }

394: /*@
395:    VecStrideNormAll - Computes the norms of subvectors of a vector defined
396:    by a starting point and a stride.

398:    Collective on Vec

400:    Input Parameter:
401: +  v - the vector
402: -  ntype - type of norm, one of NORM_1, NORM_2, NORM_INFINITY

404:    Output Parameter:
405: .  nrm - the norms

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

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

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

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

418:    Level: advanced

420:    Concepts: norm^on stride of vector
421:    Concepts: stride^norm

423: .seealso: VecNorm(), VecStrideGather(), VecStrideScatter(), VecStrideMin(), VecStrideMax()
424: @*/
425: PetscErrorCode  VecStrideNormAll(Vec v,NormType ntype,PetscReal nrm[])
426: {
427:   PetscErrorCode    ierr;
428:   PetscInt          i,j,n,bs;
429:   const PetscScalar *x;
430:   PetscReal         tnorm[128];
431:   MPI_Comm          comm;

436:   VecGetLocalSize(v,&n);
437:   VecGetArrayRead(v,&x);
438:   PetscObjectGetComm((PetscObject)v,&comm);

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

443:   if (ntype == NORM_2) {
444:     PetscScalar sum[128];
445:     for (j=0; j<bs; j++) sum[j] = 0.0;
446:     for (i=0; i<n; i+=bs) {
447:       for (j=0; j<bs; j++) sum[j] += x[i+j]*(PetscConj(x[i+j]));
448:     }
449:     for (j=0; j<bs; j++) tnorm[j]  = PetscRealPart(sum[j]);

451:     MPIU_Allreduce(tnorm,nrm,bs,MPIU_REAL,MPIU_SUM,comm);
452:     for (j=0; j<bs; j++) nrm[j] = PetscSqrtReal(nrm[j]);
453:   } else if (ntype == NORM_1) {
454:     for (j=0; j<bs; j++) tnorm[j] = 0.0;

456:     for (i=0; i<n; i+=bs) {
457:       for (j=0; j<bs; j++) tnorm[j] += PetscAbsScalar(x[i+j]);
458:     }

460:     MPIU_Allreduce(tnorm,nrm,bs,MPIU_REAL,MPIU_SUM,comm);
461:   } else if (ntype == NORM_INFINITY) {
462:     PetscReal tmp;
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++) {
467:         if ((tmp = PetscAbsScalar(x[i+j])) > tnorm[j]) tnorm[j] = tmp;
468:         /* check special case of tmp == NaN */
469:         if (tmp != tmp) {tnorm[j] = tmp; break;}
470:       }
471:     }
472:     MPIU_Allreduce(tnorm,nrm,bs,MPIU_REAL,MPIU_MAX,comm);
473:   } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown norm type");
474:   VecRestoreArrayRead(v,&x);
475:   return(0);
476: }

478: /*@
479:    VecStrideMaxAll - Computes the maximums of subvectors of a vector defined
480:    by a starting point and a stride and optionally its location.

482:    Collective on Vec

484:    Input Parameter:
485: .  v - the vector

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

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

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

498:    Level: advanced

500:    Concepts: maximum^on stride of vector
501:    Concepts: stride^maximum

503: .seealso: VecMax(), VecStrideNorm(), VecStrideGather(), VecStrideScatter(), VecStrideMin()
504: @*/
505: PetscErrorCode  VecStrideMaxAll(Vec v,PetscInt idex[],PetscReal nrm[])
506: {
507:   PetscErrorCode    ierr;
508:   PetscInt          i,j,n,bs;
509:   const PetscScalar *x;
510:   PetscReal         max[128],tmp;
511:   MPI_Comm          comm;

516:   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");
517:   VecGetLocalSize(v,&n);
518:   VecGetArrayRead(v,&x);
519:   PetscObjectGetComm((PetscObject)v,&comm);

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

524:   if (!n) {
525:     for (j=0; j<bs; j++) max[j] = PETSC_MIN_REAL;
526:   } else {
527:     for (j=0; j<bs; j++) max[j] = PetscRealPart(x[j]);

529:     for (i=bs; i<n; i+=bs) {
530:       for (j=0; j<bs; j++) {
531:         if ((tmp = PetscRealPart(x[i+j])) > max[j]) max[j] = tmp;
532:       }
533:     }
534:   }
535:   MPIU_Allreduce(max,nrm,bs,MPIU_REAL,MPIU_MAX,comm);

537:   VecRestoreArrayRead(v,&x);
538:   return(0);
539: }

541: /*@
542:    VecStrideMinAll - Computes the minimum of subvector of a vector defined
543:    by a starting point and a stride and optionally its location.

545:    Collective on Vec

547:    Input Parameter:
548: .  v - the vector

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

555:    Level: advanced

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

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

563:    Concepts: minimum^on stride of vector
564:    Concepts: stride^minimum

566: .seealso: VecMin(), VecStrideNorm(), VecStrideGather(), VecStrideScatter(), VecStrideMax()
567: @*/
568: PetscErrorCode  VecStrideMinAll(Vec v,PetscInt idex[],PetscReal nrm[])
569: {
570:   PetscErrorCode    ierr;
571:   PetscInt          i,n,bs,j;
572:   const PetscScalar *x;
573:   PetscReal         min[128],tmp;
574:   MPI_Comm          comm;

579:   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");
580:   VecGetLocalSize(v,&n);
581:   VecGetArrayRead(v,&x);
582:   PetscObjectGetComm((PetscObject)v,&comm);

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

587:   if (!n) {
588:     for (j=0; j<bs; j++) min[j] = PETSC_MAX_REAL;
589:   } else {
590:     for (j=0; j<bs; j++) min[j] = PetscRealPart(x[j]);

592:     for (i=bs; i<n; i+=bs) {
593:       for (j=0; j<bs; j++) {
594:         if ((tmp = PetscRealPart(x[i+j])) < min[j]) min[j] = tmp;
595:       }
596:     }
597:   }
598:   MPIU_Allreduce(min,nrm,bs,MPIU_REAL,MPIU_MIN,comm);

600:   VecRestoreArrayRead(v,&x);
601:   return(0);
602: }

604: /*----------------------------------------------------------------------------------------------*/
605: /*@
606:    VecStrideGatherAll - Gathers all the single components from a multi-component vector into
607:    separate vectors.

609:    Collective on Vec

611:    Input Parameter:
612: +  v - the vector
613: -  addv - one of ADD_VALUES,INSERT_VALUES,MAX_VALUES

615:    Output Parameter:
616: .  s - the location where the subvectors are stored

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

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

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

629:    Not optimized; could be easily

631:    Level: advanced

633:    Concepts: gather^into strided vector

635: .seealso: VecStrideNorm(), VecStrideScatter(), VecStrideMin(), VecStrideMax(), VecStrideGather(),
636:           VecStrideScatterAll()
637: @*/
638: PetscErrorCode  VecStrideGatherAll(Vec v,Vec s[],InsertMode addv)
639: {
640:   PetscErrorCode    ierr;
641:   PetscInt          i,n,n2,bs,j,k,*bss = NULL,nv,jj,nvc;
642:   PetscScalar       **y;
643:   const PetscScalar *x;

649:   VecGetLocalSize(v,&n);
650:   VecGetLocalSize(s[0],&n2);
651:   VecGetArrayRead(v,&x);
652:   VecGetBlockSize(v,&bs);
653:   if (bs <= 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Input vector does not have a valid blocksize set");

655:   PetscMalloc2(bs,&y,bs,&bss);
656:   nv   = 0;
657:   nvc  = 0;
658:   for (i=0; i<bs; i++) {
659:     VecGetBlockSize(s[i],&bss[i]);
660:     if (bss[i] < 1) bss[i] = 1; /* if user never set it then assume 1  Re: [PETSC #8241] VecStrideGatherAll */
661:     VecGetArray(s[i],&y[i]);
662:     nvc += bss[i];
663:     nv++;
664:     if (nvc > bs)  SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Number of subvectors in subvectors > number of vectors in main vector");
665:     if (nvc == bs) break;
666:   }

668:   n =  n/bs;

670:   jj = 0;
671:   if (addv == INSERT_VALUES) {
672:     for (j=0; j<nv; j++) {
673:       for (k=0; k<bss[j]; k++) {
674:         for (i=0; i<n; i++) y[j][i*bss[j] + k] = x[bs*i+jj+k];
675:       }
676:       jj += bss[j];
677:     }
678:   } else if (addv == ADD_VALUES) {
679:     for (j=0; j<nv; j++) {
680:       for (k=0; k<bss[j]; k++) {
681:         for (i=0; i<n; i++) y[j][i*bss[j] + k] += x[bs*i+jj+k];
682:       }
683:       jj += bss[j];
684:     }
685: #if !defined(PETSC_USE_COMPLEX)
686:   } else if (addv == MAX_VALUES) {
687:     for (j=0; j<nv; j++) {
688:       for (k=0; k<bss[j]; k++) {
689:         for (i=0; i<n; i++) y[j][i*bss[j] + k] = PetscMax(y[j][i*bss[j] + k],x[bs*i+jj+k]);
690:       }
691:       jj += bss[j];
692:     }
693: #endif
694:   } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown insert type");

696:   VecRestoreArrayRead(v,&x);
697:   for (i=0; i<nv; i++) {
698:     VecRestoreArray(s[i],&y[i]);
699:   }

701:   PetscFree2(y,bss);
702:   return(0);
703: }

705: /*@
706:    VecStrideScatterAll - Scatters all the single components from separate vectors into
707:      a multi-component vector.

709:    Collective on Vec

711:    Input Parameter:
712: +  s - the location where the subvectors are stored
713: -  addv - one of ADD_VALUES,INSERT_VALUES,MAX_VALUES

715:    Output Parameter:
716: .  v - the multicomponent vector

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

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

725:    Not optimized; could be easily

727:    Level: advanced

729:    Concepts:  scatter^into strided vector

731: .seealso: VecStrideNorm(), VecStrideScatter(), VecStrideMin(), VecStrideMax(), VecStrideGather(),
732:           VecStrideScatterAll()
733: @*/
734: PetscErrorCode  VecStrideScatterAll(Vec s[],Vec v,InsertMode addv)
735: {
736:   PetscErrorCode    ierr;
737:   PetscInt          i,n,n2,bs,j,jj,k,*bss = NULL,nv,nvc;
738:   PetscScalar       *x,**y;

744:   VecGetLocalSize(v,&n);
745:   VecGetLocalSize(s[0],&n2);
746:   VecGetArray(v,&x);
747:   VecGetBlockSize(v,&bs);
748:   if (bs <= 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Input vector does not have a valid blocksize set");

750:   PetscMalloc2(bs,&y,bs,&bss);
751:   nv   = 0;
752:   nvc  = 0;
753:   for (i=0; i<bs; i++) {
754:     VecGetBlockSize(s[i],&bss[i]);
755:     if (bss[i] < 1) bss[i] = 1; /* if user never set it then assume 1  Re: [PETSC #8241] VecStrideGatherAll */
756:     VecGetArray(s[i],&y[i]);
757:     nvc += bss[i];
758:     nv++;
759:     if (nvc > bs) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Number of subvectors in subvectors > number of vectors in main vector");
760:     if (nvc == bs) break;
761:   }

763:   n =  n/bs;

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

791:   VecRestoreArray(v,&x);
792:   for (i=0; i<nv; i++) {
793:     VecRestoreArray(s[i],&y[i]);
794:   }
795:   PetscFree2(y,bss);
796:   return(0);
797: }

799: /*@
800:    VecStrideGather - Gathers a single component from a multi-component vector into
801:    another vector.

803:    Collective on Vec

805:    Input Parameter:
806: +  v - the vector
807: .  start - starting point of the subvector (defined by a stride)
808: -  addv - one of ADD_VALUES,INSERT_VALUES,MAX_VALUES

810:    Output Parameter:
811: .  s - the location where the subvector is stored

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

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

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

823:    Not optimized; could be easily

825:    Level: advanced

827:    Concepts: gather^into strided vector

829: .seealso: VecStrideNorm(), VecStrideScatter(), VecStrideMin(), VecStrideMax(), VecStrideGatherAll(),
830:           VecStrideScatterAll()
831: @*/
832: PetscErrorCode  VecStrideGather(Vec v,PetscInt start,Vec s,InsertMode addv)
833: {

839:   if (start < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative start %D",start);
840:   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);
841:   if (!v->ops->stridegather) SETERRQ(PetscObjectComm((PetscObject)s),PETSC_ERR_SUP,"Not implemented for this Vec class");
842:   (*v->ops->stridegather)(v,start,s,addv);
843:   return(0);
844: }

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

849:    Collective on Vec

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

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

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

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

866:    Not optimized; could be easily

868:    Level: advanced

870:    Concepts: scatter^into strided vector

872: .seealso: VecStrideNorm(), VecStrideGather(), VecStrideMin(), VecStrideMax(), VecStrideGatherAll(),
873:           VecStrideScatterAll(), VecStrideSubSetScatter(), VecStrideSubSetGather()
874: @*/
875: PetscErrorCode  VecStrideScatter(Vec s,PetscInt start,Vec v,InsertMode addv)
876: {

882:   if (start < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative start %D",start);
883:   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);
884:   if (!v->ops->stridescatter) SETERRQ(PetscObjectComm((PetscObject)s),PETSC_ERR_SUP,"Not implemented for this Vec class");
885:   (*v->ops->stridescatter)(s,start,v,addv);
886:   return(0);
887: }

889: /*@
890:    VecStrideSubSetGather - Gathers a subset of components from a multi-component vector into
891:    another vector.

893:    Collective on Vec

895:    Input Parameter:
896: +  v - the vector
897: .  nidx - the number of indices
898: .  idxv - the indices of the components 0 <= idxv[0] ...idxv[nidx-1] < bs(v), they need not be sorted
899: .  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
900: -  addv - one of ADD_VALUES,INSERT_VALUES,MAX_VALUES

902:    Output Parameter:
903: .  s - the location where the subvector is stored

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


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

912:    Not optimized; could be easily

914:    Level: advanced

916:    Concepts: gather^into strided vector

918: .seealso: VecStrideNorm(), VecStrideScatter(), VecStrideGather(), VecStrideSubSetScatter(), VecStrideMin(), VecStrideMax(), VecStrideGatherAll(),
919:           VecStrideScatterAll()
920: @*/
921: PetscErrorCode  VecStrideSubSetGather(Vec v,PetscInt nidx,const PetscInt idxv[],const PetscInt idxs[],Vec s,InsertMode addv)
922: {

928:   if (nidx == PETSC_DETERMINE) nidx = s->map->bs;
929:   if (!v->ops->stridesubsetgather) SETERRQ(PetscObjectComm((PetscObject)s),PETSC_ERR_SUP,"Not implemented for this Vec class");
930:   (*v->ops->stridesubsetgather)(v,nidx,idxv,idxs,s,addv);
931:   return(0);
932: }

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

937:    Collective on Vec

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

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

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

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

955:    Not optimized; could be easily

957:    Level: advanced

959:    Concepts: scatter^into strided vector

961: .seealso: VecStrideNorm(), VecStrideGather(), VecStrideGather(), VecStrideSubSetGather(), VecStrideMin(), VecStrideMax(), VecStrideGatherAll(),
962:           VecStrideScatterAll()
963: @*/
964: PetscErrorCode  VecStrideSubSetScatter(Vec s,PetscInt nidx,const PetscInt idxs[],const PetscInt idxv[],Vec v,InsertMode addv)
965: {

971:   if (nidx == PETSC_DETERMINE) nidx = s->map->bs;
972:   if (!v->ops->stridesubsetscatter) SETERRQ(PetscObjectComm((PetscObject)s),PETSC_ERR_SUP,"Not implemented for this Vec class");
973:   (*v->ops->stridesubsetscatter)(s,nidx,idxs,idxv,v,addv);
974:   return(0);
975: }

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

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

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

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

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

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

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

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

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

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

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

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

1056:   bs  = v->map->bs;
1057:   bss = s->map->bs;
1058:   n  =  n/bs;

1060: #if defined(PETSC_DEBUG)
1061:   if (n != ns/bss) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Incompatible layout of vectors");
1062:   for (j=0; j<nidx; j++) {
1063:     if (idxv[j] < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"idx[%D] %D is negative",j,idxv[j]);
1064:     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);
1065:   }
1066:   if (!idxs && bss != nidx) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Must provide idxs when not gathering into all locations");
1067: #endif

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

1103:   VecRestoreArrayRead(v,&x);
1104:   VecRestoreArray(s,&y);
1105:   return(0);
1106: }

1108: PetscErrorCode  VecStrideSubSetScatter_Default(Vec s,PetscInt nidx,const PetscInt idxs[],const PetscInt idxv[],Vec v,InsertMode addv)
1109: {
1110:   PetscErrorCode    ierr;
1111:   PetscInt          j,i,n,bs,ns,bss;
1112:   PetscScalar       *x;
1113:   const PetscScalar *y;

1116:   VecGetLocalSize(v,&n);
1117:   VecGetLocalSize(s,&ns);
1118:   VecGetArray(v,&x);
1119:   VecGetArrayRead(s,&y);

1121:   bs  = v->map->bs;
1122:   bss = s->map->bs;
1123:   n  =  n/bs;

1125: #if defined(PETSC_DEBUG)
1126:   if (n != ns/bss) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Incompatible layout of vectors");
1127:   for (j=0; j<bss; j++) {
1128:     if (idx[j] < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"idx[%D] %D is negative",j,idx[j]);
1129:     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);
1130:   }
1131:   if (!idxs && bss != nidx) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Must provide idxs when not scattering from all locations");
1132: #endif

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

1168:   VecRestoreArray(v,&x);
1169:   VecRestoreArrayRead(s,&y);
1170:   return(0);
1171: }

1173: PetscErrorCode VecReciprocal_Default(Vec v)
1174: {
1176:   PetscInt       i,n;
1177:   PetscScalar    *x;

1180:   VecGetLocalSize(v,&n);
1181:   VecGetArray(v,&x);
1182:   for (i=0; i<n; i++) {
1183:     if (x[i] != (PetscScalar)0.0) x[i] = (PetscScalar)1.0/x[i];
1184:   }
1185:   VecRestoreArray(v,&x);
1186:   return(0);
1187: }

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

1192:   Not collective

1194:   Input Parameter:
1195: . v - The vector

1197:   Output Parameter:
1198: . v - The vector of exponents

1200:   Level: beginner

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

1204: .keywords: vector, sqrt, square root
1205: @*/
1206: PetscErrorCode  VecExp(Vec v)
1207: {
1208:   PetscScalar    *x;
1209:   PetscInt       i, n;

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

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

1228:   Not collective

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

1233:   Output Parameter:
1234: . v - The vector of logs

1236:   Level: beginner

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

1240: .keywords: vector, sqrt, square root
1241: @*/
1242: PetscErrorCode  VecLog(Vec v)
1243: {
1244:   PetscScalar    *x;
1245:   PetscInt       i, n;

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

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

1264:   Not collective

1266:   Input Parameter:
1267: . v - The vector

1269:   Output Parameter:
1270: . v - The vector square root

1272:   Level: beginner

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

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

1278: .keywords: vector, sqrt, square root
1279: @*/
1280: PetscErrorCode  VecSqrtAbs(Vec v)
1281: {
1282:   PetscScalar    *x;
1283:   PetscInt       i, n;

1288:   if (v->ops->sqrt) {
1289:     (*v->ops->sqrt)(v);
1290:   } else {
1291:     VecGetLocalSize(v, &n);
1292:     VecGetArray(v, &x);
1293:     for (i = 0; i < n; i++) x[i] = PetscSqrtReal(PetscAbsScalar(x[i]));
1294:     VecRestoreArray(v, &x);
1295:   }
1296:   return(0);
1297: }

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

1302:   Collective on Vec

1304:   Input Parameter:
1305: + s - first vector
1306: - t - second vector

1308:   Output Parameter:
1309: + dp - s'conj(t)
1310: - nm - t'conj(t)

1312:   Level: advanced

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


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

1319: .keywords: vector, sqrt, square root
1320: @*/
1321: PetscErrorCode  VecDotNorm2(Vec s,Vec t,PetscScalar *dp, PetscReal *nm)
1322: {
1323:   const PetscScalar *sx, *tx;
1324:   PetscScalar       dpx = 0.0, nmx = 0.0,work[2],sum[2];
1325:   PetscInt          i, n;
1326:   PetscErrorCode    ierr;

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

1339:   PetscLogEventBarrierBegin(VEC_DotNormBarrier,s,t,0,0,PetscObjectComm((PetscObject)s));
1340:   if (s->ops->dotnorm2) {
1341:     (*s->ops->dotnorm2)(s,t,dp,&dpx);
1342:     *nm  = PetscRealPart(dpx);
1343:   } else {
1344:     VecGetLocalSize(s, &n);
1345:     VecGetArrayRead(s, &sx);
1346:     VecGetArrayRead(t, &tx);

1348:     for (i = 0; i<n; i++) {
1349:       dpx += sx[i]*PetscConj(tx[i]);
1350:       nmx += tx[i]*PetscConj(tx[i]);
1351:     }
1352:     work[0] = dpx;
1353:     work[1] = nmx;

1355:     MPIU_Allreduce(work,sum,2,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)s));
1356:     *dp  = sum[0];
1357:     *nm  = PetscRealPart(sum[1]);

1359:     VecRestoreArrayRead(t, &tx);
1360:     VecRestoreArrayRead(s, &sx);
1361:     PetscLogFlops(4.0*n);
1362:   }
1363:   PetscLogEventBarrierEnd(VEC_DotNormBarrier,s,t,0,0,PetscObjectComm((PetscObject)s));
1364:   return(0);
1365: }

1367: /*@
1368:    VecSum - Computes the sum of all the components of a vector.

1370:    Collective on Vec

1372:    Input Parameter:
1373: .  v - the vector

1375:    Output Parameter:
1376: .  sum - the result

1378:    Level: beginner

1380:    Concepts: sum^of vector entries

1382: .seealso: VecNorm()
1383: @*/
1384: PetscErrorCode  VecSum(Vec v,PetscScalar *sum)
1385: {
1386:   PetscErrorCode    ierr;
1387:   PetscInt          i,n;
1388:   const PetscScalar *x;
1389:   PetscScalar       lsum = 0.0;

1394:   VecGetLocalSize(v,&n);
1395:   VecGetArrayRead(v,&x);
1396:   for (i=0; i<n; i++) lsum += x[i];
1397:   MPIU_Allreduce(&lsum,sum,1,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)v));
1398:   VecRestoreArrayRead(v,&x);
1399:   return(0);
1400: }

1402: /*@
1403:    VecShift - Shifts all of the components of a vector by computing
1404:    x[i] = x[i] + shift.

1406:    Logically Collective on Vec

1408:    Input Parameters:
1409: +  v - the vector
1410: -  shift - the shift

1412:    Output Parameter:
1413: .  v - the shifted vector

1415:    Level: intermediate

1417:    Concepts: vector^adding constant

1419: @*/
1420: PetscErrorCode  VecShift(Vec v,PetscScalar shift)
1421: {
1423:   PetscInt       i,n;
1424:   PetscScalar    *x;

1429:   VecLocked(v,1);

1431:   if (v->ops->shift) {
1432:     (*v->ops->shift)(v,shift);
1433:   } else {
1434:     VecGetLocalSize(v,&n);
1435:     VecGetArray(v,&x);
1436:     for (i=0; i<n; i++) x[i] += shift;
1437:     VecRestoreArray(v,&x);
1438:   }
1439:   return(0);
1440: }

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

1445:    Logically Collective on Vec

1447:    Input Parameters:
1448: .  v - the vector

1450:    Level: intermediate

1452:    Concepts: vector^absolute value

1454: @*/
1455: PetscErrorCode  VecAbs(Vec v)
1456: {
1458:   PetscInt       i,n;
1459:   PetscScalar    *x;

1463:   VecLocked(v,1);
1464: 
1465:   if (v->ops->abs) {
1466:     (*v->ops->abs)(v);
1467:   } else {
1468:     VecGetLocalSize(v,&n);
1469:     VecGetArray(v,&x);
1470:     for (i=0; i<n; i++) x[i] = PetscAbsScalar(x[i]);
1471:     VecRestoreArray(v,&x);
1472:   }
1473:   return(0);
1474: }

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

1479:   Input Parameters:
1480: + vec   - The vector
1481: . order - The ordering
1482: - inv   - The flag for inverting the permutation

1484:   Level: beginner

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

1488: .seealso: MatPermute()
1489: .keywords: vec, permute
1490: @*/
1491: PetscErrorCode  VecPermute(Vec x, IS row, PetscBool inv)
1492: {
1493:   PetscScalar    *array, *newArray;
1494:   const PetscInt *idx;
1495:   PetscInt       i,rstart,rend;

1499:   VecLocked(x,1);
1500: 
1501:   VecGetOwnershipRange(x,&rstart,&rend);
1502:   ISGetIndices(row, &idx);
1503:   VecGetArray(x, &array);
1504:   PetscMalloc1(x->map->n, &newArray);
1505: #if defined(PETSC_USE_DEBUG)
1506:   for (i = 0; i < x->map->n; i++) {
1507:     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]);
1508:   }
1509: #endif
1510:   if (!inv) {
1511:     for (i = 0; i < x->map->n; i++) newArray[i] = array[idx[i]-rstart];
1512:   } else {
1513:     for (i = 0; i < x->map->n; i++) newArray[idx[i]-rstart] = array[i];
1514:   }
1515:   VecRestoreArray(x, &array);
1516:   ISRestoreIndices(row, &idx);
1517:   VecReplaceArray(x, newArray);
1518:   return(0);
1519: }

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

1526:    Collective on Vec

1528:    Input Parameters:
1529: +  vec1 - the first vector
1530: -  vec2 - the second vector

1532:    Output Parameter:
1533: .  flg - PETSC_TRUE if the vectors are equal; PETSC_FALSE otherwise.

1535:    Level: intermediate

1537:    Concepts: equal^two vectors
1538:    Concepts: vector^equality

1540: @*/
1541: PetscErrorCode  VecEqual(Vec vec1,Vec vec2,PetscBool  *flg)
1542: {
1543:   const PetscScalar  *v1,*v2;
1544:   PetscErrorCode     ierr;
1545:   PetscInt           n1,n2,N1,N2;
1546:   PetscBool          flg1;

1552:   if (vec1 == vec2) *flg = PETSC_TRUE;
1553:   else {
1554:     VecGetSize(vec1,&N1);
1555:     VecGetSize(vec2,&N2);
1556:     if (N1 != N2) flg1 = PETSC_FALSE;
1557:     else {
1558:       VecGetLocalSize(vec1,&n1);
1559:       VecGetLocalSize(vec2,&n2);
1560:       if (n1 != n2) flg1 = PETSC_FALSE;
1561:       else {
1562:         VecGetArrayRead(vec1,&v1);
1563:         VecGetArrayRead(vec2,&v2);
1564:         PetscMemcmp(v1,v2,n1*sizeof(PetscScalar),&flg1);
1565:         VecRestoreArrayRead(vec1,&v1);
1566:         VecRestoreArrayRead(vec2,&v2);
1567:       }
1568:     }
1569:     /* combine results from all processors */
1570:     MPIU_Allreduce(&flg1,flg,1,MPIU_BOOL,MPI_MIN,PetscObjectComm((PetscObject)vec1));
1571:   }
1572:   return(0);
1573: }

1575: /*@
1576:    VecUniqueEntries - Compute the number of unique entries, and those entries

1578:    Collective on Vec

1580:    Input Parameter:
1581: .  vec - the vector

1583:    Output Parameters:
1584: +  n - The number of unique entries
1585: -  e - The entries

1587:    Level: intermediate

1589: @*/
1590: PetscErrorCode  VecUniqueEntries(Vec vec, PetscInt *n, PetscScalar **e)
1591: {
1592:   const PetscScalar *v;
1593:   PetscScalar       *tmp, *vals;
1594:   PetscMPIInt       *N, *displs, l;
1595:   PetscInt          ng, m, i, j, p;
1596:   PetscMPIInt       size;
1597:   PetscErrorCode    ierr;

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