Actual source code: vinv.c

petsc-master 2019-10-20
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


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

 41:   VecSetErrorIfLocked(v,1);

 43:   VecGetLocalSize(v,&n);
 44:   VecGetArray(v,&x);
 45:   VecGetBlockSize(v,&bs);
 46:   if (start < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative start %D",start);
 47:   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);
 48:   x += start;
 49:   for (i=0; i<n; i+=bs) x[i] = s;
 50:   x -= start;
 51:   VecRestoreArray(v,&x);
 52:   return(0);
 53: }

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

 59:    Logically Collective on Vec

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

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

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

 72:    Level: advanced


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

 87:   VecSetErrorIfLocked(v,1);

 89:   VecGetLocalSize(v,&n);
 90:   VecGetArray(v,&x);
 91:   VecGetBlockSize(v,&bs);
 92:   if (start < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative start %D",start);
 93:   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);
 94:   x += start;
 95:   for (i=0; i<n; i+=bs) x[i] *= scale;
 96:   x -= start;
 97:   VecRestoreArray(v,&x);
 98:   return(0);
 99: }

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

105:    Collective on Vec

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

112:    Output Parameter:
113: .  norm - the norm

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

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

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

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

127:    Level: advanced


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

143:   VecGetLocalSize(v,&n);
144:   VecGetArrayRead(v,&x);
145:   PetscObjectGetComm((PetscObject)v,&comm);

147:   VecGetBlockSize(v,&bs);
148:   if (start < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative start %D",start);
149:   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);
150:   x += start;

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

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

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

181:    Collective on Vec

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

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

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

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

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

202:    Level: advanced


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


219:   VecGetLocalSize(v,&n);
220:   VecGetArrayRead(v,&x);
221:   PetscObjectGetComm((PetscObject)v,&comm);

223:   VecGetBlockSize(v,&bs);
224:   if (start < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative start %D",start);
225:   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);
226:   x += start;

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

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

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

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

259:    Collective on Vec

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

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

269:    Level: advanced

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

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

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


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


297:   VecGetLocalSize(v,&n);
298:   VecGetArrayRead(v,&x);
299:   PetscObjectGetComm((PetscObject)v,&comm);

301:   VecGetBlockSize(v,&bs);
302:   if (start < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative start %D",start);
303:   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);
304:   x += start;

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

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

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

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

337:    Logically Collective on Vec

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

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

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


350:    Level: advanced


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

364:   VecSetErrorIfLocked(v,1);

366:   VecGetLocalSize(v,&n);
367:   VecGetArray(v,&x);
368:   VecGetBlockSize(v,&bs);

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

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

382:    Collective on Vec

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

388:    Output Parameter:
389: .  nrm - the norms

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

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

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

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

402:    Level: advanced


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

418:   VecGetLocalSize(v,&n);
419:   VecGetArrayRead(v,&x);
420:   PetscObjectGetComm((PetscObject)v,&comm);

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

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

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

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

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

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

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

464:    Collective on Vec

466:    Input Parameter:
467: .  v - the vector

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

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

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

480:    Level: advanced


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

496:   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");
497:   VecGetLocalSize(v,&n);
498:   VecGetArrayRead(v,&x);
499:   PetscObjectGetComm((PetscObject)v,&comm);

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

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

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

517:   VecRestoreArrayRead(v,&x);
518:   return(0);
519: }

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

525:    Collective on Vec

527:    Input Parameter:
528: .  v - the vector

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

535:    Level: advanced

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

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


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

557:   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");
558:   VecGetLocalSize(v,&n);
559:   VecGetArrayRead(v,&x);
560:   PetscObjectGetComm((PetscObject)v,&comm);

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

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

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

578:   VecRestoreArrayRead(v,&x);
579:   return(0);
580: }

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

587:    Collective on Vec

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

593:    Output Parameter:
594: .  s - the location where the subvectors are stored

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

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

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

607:    Not optimized; could be easily

609:    Level: advanced

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

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

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

644:   n =  n/bs;

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

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

677:   PetscFree2(y,bss);
678:   return(0);
679: }

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

685:    Collective on Vec

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

691:    Output Parameter:
692: .  v - the multicomponent vector

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

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

701:    Not optimized; could be easily

703:    Level: advanced

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

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

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

738:   n =  n/bs;

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

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

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

778:    Collective on Vec

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

785:    Output Parameter:
786: .  s - the location where the subvector is stored

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

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

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

798:    Not optimized; could be easily

800:    Level: advanced

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

812:   if (start < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative start %D",start);
813:   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);
814:   if (!v->ops->stridegather) SETERRQ(PetscObjectComm((PetscObject)s),PETSC_ERR_SUP,"Not implemented for this Vec class");
815:   (*v->ops->stridegather)(v,start,s,addv);
816:   return(0);
817: }

819: /*@
820:    VecStrideScatter - Scatters a single component from a vector into a multi-component vector.

822:    Collective on Vec

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

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

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

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

839:    Not optimized; could be easily

841:    Level: advanced

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

853:   if (start < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative start %D",start);
854:   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);
855:   if (!v->ops->stridescatter) SETERRQ(PetscObjectComm((PetscObject)s),PETSC_ERR_SUP,"Not implemented for this Vec class");
856:   (*v->ops->stridescatter)(s,start,v,addv);
857:   return(0);
858: }

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

864:    Collective on Vec

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

873:    Output Parameter:
874: .  s - the location where the subvector is stored

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


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

883:    Not optimized; could be easily

885:    Level: advanced

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

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

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

906:    Collective on Vec

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

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

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

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

924:    Not optimized; could be easily

926:    Level: advanced

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

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

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

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

957:   bs = v->map->bs;
958:   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);
959:   x += start;
960:   n  =  n/bs;

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

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

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

985:   VecGetLocalSize(v,&n);
986:   VecGetLocalSize(s,&ns);
987:   VecGetArray(v,&x);
988:   VecGetArrayRead(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 scatter to multicomponent vector %D",ns*bs,n);
992:   x += start;
993:   n  =  n/bs;

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

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

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

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

1023:   bs  = v->map->bs;
1024:   bss = s->map->bs;
1025:   n  =  n/bs;

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

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

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

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

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

1088:   bs  = v->map->bs;
1089:   bss = s->map->bs;
1090:   n  =  n/bs;

1092: #if defined(PETSC_USE_DEBUG)
1093:   if (n != ns/bss) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Incompatible layout of vectors");
1094:   for (j=0; j<bss; j++) {
1095:     if (idxs) {
1096:       if (idxs[j] < 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"idx[%D] %D is negative",j,idxs[j]);
1097:       if (idxs[j] >= bs) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"idx[%D] %D is greater than or equal to vector blocksize %D",j,idxs[j],bs);
1098:     }
1099:   }
1100:   if (!idxs && bss != nidx) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Must provide idxs when not scattering from 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++) x[bs*i + idxv[j]] = y[bss*i+j];
1107:       }
1108:     } else {
1109:       for (i=0; i<n; i++) {
1110:         for (j=0; j<bss; j++) x[bs*i + idxv[j]] = y[bss*i+idxs[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++) x[bs*i + idxv[j]] += y[bss*i+j];
1117:       }
1118:     } else {
1119:       for (i=0; i<n; i++) {
1120:         for (j=0; j<bss; j++) x[bs*i + idxv[j]] += y[bss*i+idxs[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++) x[bs*i + idxv[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++) x[bs*i + idxv[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:   VecRestoreArray(v,&x);
1138:   VecRestoreArrayRead(s,&y);
1139:   return(0);
1140: }

1142: PetscErrorCode VecReciprocal_Default(Vec v)
1143: {
1145:   PetscInt       i,n;
1146:   PetscScalar    *x;

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

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

1161:   Not collective

1163:   Input Parameter:
1164: . v - The vector

1166:   Output Parameter:
1167: . v - The vector of exponents

1169:   Level: beginner

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

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

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

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

1196:   Not collective

1198:   Input Parameter:
1199: . v - The vector

1201:   Output Parameter:
1202: . v - The vector of logs

1204:   Level: beginner

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

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

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

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

1231:   Not collective

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

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

1239:   Level: beginner

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

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

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

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

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

1268:   Collective on Vec

1270:   Input Parameter:
1271: + s - first vector
1272: - t - second vector

1274:   Output Parameter:
1275: + dp - s'conj(t)
1276: - nm - t'conj(t)

1278:   Level: advanced

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


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

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

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

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

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

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

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

1333: /*@
1334:    VecSum - Computes the sum of all the components of a vector.

1336:    Collective on Vec

1338:    Input Parameter:
1339: .  v - the vector

1341:    Output Parameter:
1342: .  sum - the result

1344:    Level: beginner

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

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

1366: /*@
1367:    VecImaginaryPart - Replaces a complex vector with its imginary part

1369:    Collective on Vec

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

1374:    Level: beginner

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

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

1393: /*@
1394:    VecRealPart - Replaces a complex vector with its real part

1396:    Collective on Vec

1398:    Input Parameter:
1399: .  v - the vector

1401:    Level: beginner

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

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

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

1424:    Logically Collective on Vec

1426:    Input Parameters:
1427: +  v - the vector
1428: -  shift - the shift

1430:    Output Parameter:
1431: .  v - the shifted vector

1433:    Level: intermediate

1435: @*/
1436: PetscErrorCode  VecShift(Vec v,PetscScalar shift)
1437: {
1439:   PetscInt       i,n;
1440:   PetscScalar    *x;

1445:   VecSetErrorIfLocked(v,1);
1446:   if (shift == 0.0) return(0);

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

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

1462:    Logically Collective on Vec

1464:    Input Parameters:
1465: .  v - the vector

1467:    Level: intermediate

1469: @*/
1470: PetscErrorCode  VecAbs(Vec v)
1471: {
1473:   PetscInt       i,n;
1474:   PetscScalar    *x;

1478:   VecSetErrorIfLocked(v,1);

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

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

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

1499:   Level: beginner

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

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

1514:   VecSetErrorIfLocked(x,1);

1516:   VecGetOwnershipRange(x,&rstart,&rend);
1517:   ISGetIndices(row, &idx);
1518:   VecGetArrayRead(x, &array);
1519:   PetscMalloc1(x->map->n, &newArray);
1520: #if defined(PETSC_USE_DEBUG)
1521:   for (i = 0; i < x->map->n; i++) {
1522:     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]);
1523:   }
1524: #endif
1525:   if (!inv) {
1526:     for (i = 0; i < x->map->n; i++) newArray[i] = array[idx[i]-rstart];
1527:   } else {
1528:     for (i = 0; i < x->map->n; i++) newArray[idx[i]-rstart] = array[i];
1529:   }
1530:   VecRestoreArrayRead(x, &array);
1531:   ISRestoreIndices(row, &idx);
1532:   VecReplaceArray(x, newArray);
1533:   return(0);
1534: }

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

1541:    Collective on Vec

1543:    Input Parameters:
1544: +  vec1 - the first vector
1545: -  vec2 - the second vector

1547:    Output Parameter:
1548: .  flg - PETSC_TRUE if the vectors are equal; PETSC_FALSE otherwise.

1550:    Level: intermediate


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

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

1588: /*@
1589:    VecUniqueEntries - Compute the number of unique entries, and those entries

1591:    Collective on Vec

1593:    Input Parameter:
1594: .  vec - the vector

1596:    Output Parameters:
1597: +  n - The number of unique entries
1598: -  e - The entries

1600:    Level: intermediate

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

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