Actual source code: vinv.c


  2: /*
  3:      Some useful vector utility functions.
  4: */
  5: #include <../src/vec/vec/impls/mpi/pvecimpl.h>

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

 11:    Logically Collective on Vec

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

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

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

 24:    Level: advanced

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

 38:   VecSetErrorIfLocked(v,1);

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

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

 56:    Logically Collective on Vec

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

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

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

 69:    Level: advanced

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

 83:   VecSetErrorIfLocked(v,1);

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

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

101:    Collective on Vec

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

108:    Output Parameter:
109: .  norm - the norm

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

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

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

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

123:    Level: advanced

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

138:   VecGetLocalSize(v,&n);
139:   VecGetArrayRead(v,&x);
140:   PetscObjectGetComm((PetscObject)v,&comm);

142:   VecGetBlockSize(v,&bs);
143:   if (start < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative start %D",start);
144:   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);
145:   x += start;

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

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

172: /*@
173:    VecStrideMax - Computes the maximum of subvector of a vector defined
174:    by a starting point and a stride and optionally its location.

176:    Collective on Vec

178:    Input Parameters:
179: +  v - the vector
180: -  start - starting point of the subvector (defined by a stride)

182:    Output Parameters:
183: +  idex - the location where the maximum occurred  (pass NULL if not required)
184: -  nrm - the maximum value in the subvector

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

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

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

197:    Level: advanced

199: .seealso: VecMax(), VecStrideNorm(), VecStrideGather(), VecStrideScatter(), VecStrideMin()
200: @*/
201: PetscErrorCode  VecStrideMax(Vec v,PetscInt start,PetscInt *idex,PetscReal *nrm)
202: {
203:   PetscErrorCode    ierr;
204:   PetscInt          i,n,bs,id;
205:   const PetscScalar *x;
206:   PetscReal         max,tmp;


212:   VecGetLocalSize(v,&n);
213:   VecGetArrayRead(v,&x);

215:   VecGetBlockSize(v,&bs);
216:   if (start < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative start %D",start);
217:   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);
218:   x += start;

220:   id = -1;
221:   if (!n) max = PETSC_MIN_REAL;
222:   else {
223:     id  = 0;
224:     max = PetscRealPart(x[0]);
225:     for (i=bs; i<n; i+=bs) {
226:       if ((tmp = PetscRealPart(x[i])) > max) { max = tmp; id = i;}
227:     }
228:   }
229:   VecRestoreArrayRead(v,&x);

231:   if (!idex) {
232:     MPIU_Allreduce(&max,nrm,1,MPIU_REAL,MPIU_MAX,PetscObjectComm((PetscObject)v));
233:   } else {
234:     struct { PetscReal v; PetscInt i; } in,out;
235:     PetscInt rstart;

237:     VecGetOwnershipRange(v,&rstart,NULL);
238:     in.v  = max;
239:     in.i  = rstart+id+start;
240:     MPIU_Allreduce(&in,&out,1,MPIU_REAL_INT,MPIU_MAXLOC,PetscObjectComm((PetscObject)v));
241:     *nrm  = out.v;
242:     *idex = out.i;
243:   }
244:   return(0);
245: }

247: /*@
248:    VecStrideMin - Computes the minimum of subvector of a vector defined
249:    by a starting point and a stride and optionally its location.

251:    Collective on Vec

253:    Input Parameters:
254: +  v - the vector
255: -  start - starting point of the subvector (defined by a stride)

257:    Output Parameters:
258: +  idex - the location where the minimum occurred. (pass NULL if not required)
259: -  nrm - the minimum value in the subvector

261:    Level: advanced

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

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

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

274: .seealso: VecMin(), VecStrideNorm(), VecStrideGather(), VecStrideScatter(), VecStrideMax()
275: @*/
276: PetscErrorCode  VecStrideMin(Vec v,PetscInt start,PetscInt *idex,PetscReal *nrm)
277: {
278:   PetscErrorCode    ierr;
279:   PetscInt          i,n,bs,id;
280:   const PetscScalar *x;
281:   PetscReal         min,tmp;
282:   MPI_Comm          comm;


288:   VecGetLocalSize(v,&n);
289:   VecGetArrayRead(v,&x);
290:   PetscObjectGetComm((PetscObject)v,&comm);

292:   VecGetBlockSize(v,&bs);
293:   if (start < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative start %D",start);
294:   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);
295:   x += start;

297:   id = -1;
298:   if (!n) min = PETSC_MAX_REAL;
299:   else {
300:     id = 0;
301:     min = PetscRealPart(x[0]);
302:     for (i=bs; i<n; i+=bs) {
303:       if ((tmp = PetscRealPart(x[i])) < min) { min = tmp; id = i;}
304:     }
305:   }
306:   VecRestoreArrayRead(v,&x);

308:   if (!idex) {
309:     MPIU_Allreduce(&min,nrm,1,MPIU_REAL,MPIU_MIN,comm);
310:   } else {
311:     struct { PetscReal v; PetscInt i; } in,out;
312:     PetscInt rstart;

314:     VecGetOwnershipRange(v,&rstart,NULL);
315:     in.v  = min;
316:     in.i  = rstart+id+start;
317:     MPIU_Allreduce(&in,&out,1,MPIU_REAL_INT,MPIU_MINLOC,PetscObjectComm((PetscObject)v));
318:     *nrm  = out.v;
319:     *idex = out.i;
320:   }
321:   return(0);
322: }

324: /*@
325:    VecStrideScaleAll - Scales the subvectors of a vector defined
326:    by a starting point and a stride.

328:    Logically Collective on Vec

330:    Input Parameters:
331: +  v - the vector
332: -  scales - values to multiply each subvector entry by

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

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

340:    Level: advanced

342: .seealso: VecNorm(), VecStrideScale(), VecScale(), VecStrideGather(), VecStrideScatter(), VecStrideMin(), VecStrideMax()
343: @*/
344: PetscErrorCode  VecStrideScaleAll(Vec v,const PetscScalar *scales)
345: {
347:   PetscInt       i,j,n,bs;
348:   PetscScalar    *x;

353:   VecSetErrorIfLocked(v,1);

355:   VecGetLocalSize(v,&n);
356:   VecGetArray(v,&x);
357:   VecGetBlockSize(v,&bs);

359:   /* need to provide optimized code for each bs */
360:   for (i=0; i<n; i+=bs) {
361:     for (j=0; j<bs; j++) x[i+j] *= scales[j];
362:   }
363:   VecRestoreArray(v,&x);
364:   return(0);
365: }

367: /*@
368:    VecStrideNormAll - Computes the norms of subvectors of a vector defined
369:    by a starting point and a stride.

371:    Collective on Vec

373:    Input Parameters:
374: +  v - the vector
375: -  ntype - type of norm, one of NORM_1, NORM_2, NORM_INFINITY

377:    Output Parameter:
378: .  nrm - the norms

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

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

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

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

391:    Level: advanced

393: .seealso: VecNorm(), VecStrideGather(), VecStrideScatter(), VecStrideMin(), VecStrideMax()
394: @*/
395: PetscErrorCode  VecStrideNormAll(Vec v,NormType ntype,PetscReal nrm[])
396: {
397:   PetscErrorCode    ierr;
398:   PetscInt          i,j,n,bs;
399:   const PetscScalar *x;
400:   PetscReal         tnorm[128];
401:   MPI_Comm          comm;

406:   VecGetLocalSize(v,&n);
407:   VecGetArrayRead(v,&x);
408:   PetscObjectGetComm((PetscObject)v,&comm);

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

413:   if (ntype == NORM_2) {
414:     PetscScalar sum[128];
415:     for (j=0; j<bs; j++) sum[j] = 0.0;
416:     for (i=0; i<n; i+=bs) {
417:       for (j=0; j<bs; j++) sum[j] += x[i+j]*(PetscConj(x[i+j]));
418:     }
419:     for (j=0; j<bs; j++) tnorm[j]  = PetscRealPart(sum[j]);

421:     MPIU_Allreduce(tnorm,nrm,bs,MPIU_REAL,MPIU_SUM,comm);
422:     for (j=0; j<bs; j++) nrm[j] = PetscSqrtReal(nrm[j]);
423:   } else if (ntype == NORM_1) {
424:     for (j=0; j<bs; j++) tnorm[j] = 0.0;

426:     for (i=0; i<n; i+=bs) {
427:       for (j=0; j<bs; j++) tnorm[j] += PetscAbsScalar(x[i+j]);
428:     }

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

435:     for (i=0; i<n; i+=bs) {
436:       for (j=0; j<bs; j++) {
437:         if ((tmp = PetscAbsScalar(x[i+j])) > tnorm[j]) tnorm[j] = tmp;
438:         /* check special case of tmp == NaN */
439:         if (tmp != tmp) {tnorm[j] = tmp; break;}
440:       }
441:     }
442:     MPIU_Allreduce(tnorm,nrm,bs,MPIU_REAL,MPIU_MAX,comm);
443:   } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown norm type");
444:   VecRestoreArrayRead(v,&x);
445:   return(0);
446: }

448: /*@
449:    VecStrideMaxAll - Computes the maximums of subvectors of a vector defined
450:    by a starting point and a stride and optionally its location.

452:    Collective on Vec

454:    Input Parameter:
455: .  v - the vector

457:    Output Parameters:
458: +  index - the location where the maximum occurred (not supported, pass NULL,
459:            if you need this, send mail to petsc-maint@mcs.anl.gov to request it)
460: -  nrm - the maximum values of each subvector

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

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

468:    Level: advanced

470: .seealso: VecMax(), VecStrideNorm(), VecStrideGather(), VecStrideScatter(), VecStrideMin()
471: @*/
472: PetscErrorCode  VecStrideMaxAll(Vec v,PetscInt idex[],PetscReal nrm[])
473: {
474:   PetscErrorCode    ierr;
475:   PetscInt          i,j,n,bs;
476:   const PetscScalar *x;
477:   PetscReal         max[128],tmp;
478:   MPI_Comm          comm;

483:   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");
484:   VecGetLocalSize(v,&n);
485:   VecGetArrayRead(v,&x);
486:   PetscObjectGetComm((PetscObject)v,&comm);

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

491:   if (!n) {
492:     for (j=0; j<bs; j++) max[j] = PETSC_MIN_REAL;
493:   } else {
494:     for (j=0; j<bs; j++) max[j] = PetscRealPart(x[j]);

496:     for (i=bs; i<n; i+=bs) {
497:       for (j=0; j<bs; j++) {
498:         if ((tmp = PetscRealPart(x[i+j])) > max[j]) max[j] = tmp;
499:       }
500:     }
501:   }
502:   MPIU_Allreduce(max,nrm,bs,MPIU_REAL,MPIU_MAX,comm);

504:   VecRestoreArrayRead(v,&x);
505:   return(0);
506: }

508: /*@
509:    VecStrideMinAll - Computes the minimum of subvector of a vector defined
510:    by a starting point and a stride and optionally its location.

512:    Collective on Vec

514:    Input Parameter:
515: .  v - the vector

517:    Output Parameters:
518: +  idex - the location where the minimum occurred (not supported, pass NULL,
519:            if you need this, send mail to petsc-maint@mcs.anl.gov to request it)
520: -  nrm - the minimums of each subvector

522:    Level: advanced

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

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

530: .seealso: VecMin(), VecStrideNorm(), VecStrideGather(), VecStrideScatter(), VecStrideMax()
531: @*/
532: PetscErrorCode  VecStrideMinAll(Vec v,PetscInt idex[],PetscReal nrm[])
533: {
534:   PetscErrorCode    ierr;
535:   PetscInt          i,n,bs,j;
536:   const PetscScalar *x;
537:   PetscReal         min[128],tmp;
538:   MPI_Comm          comm;

543:   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");
544:   VecGetLocalSize(v,&n);
545:   VecGetArrayRead(v,&x);
546:   PetscObjectGetComm((PetscObject)v,&comm);

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

551:   if (!n) {
552:     for (j=0; j<bs; j++) min[j] = PETSC_MAX_REAL;
553:   } else {
554:     for (j=0; j<bs; j++) min[j] = PetscRealPart(x[j]);

556:     for (i=bs; i<n; i+=bs) {
557:       for (j=0; j<bs; j++) {
558:         if ((tmp = PetscRealPart(x[i+j])) < min[j]) min[j] = tmp;
559:       }
560:     }
561:   }
562:   MPIU_Allreduce(min,nrm,bs,MPIU_REAL,MPIU_MIN,comm);

564:   VecRestoreArrayRead(v,&x);
565:   return(0);
566: }

568: /*----------------------------------------------------------------------------------------------*/
569: /*@
570:    VecStrideGatherAll - Gathers all the single components from a multi-component vector into
571:    separate vectors.

573:    Collective on Vec

575:    Input Parameters:
576: +  v - the vector
577: -  addv - one of ADD_VALUES,INSERT_VALUES,MAX_VALUES

579:    Output Parameter:
580: .  s - the location where the subvectors are stored

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

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

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

593:    Not optimized; could be easily

595:    Level: advanced

597: .seealso: VecStrideNorm(), VecStrideScatter(), VecStrideMin(), VecStrideMax(), VecStrideGather(),
598:           VecStrideScatterAll()
599: @*/
600: PetscErrorCode  VecStrideGatherAll(Vec v,Vec s[],InsertMode addv)
601: {
602:   PetscErrorCode    ierr;
603:   PetscInt          i,n,n2,bs,j,k,*bss = NULL,nv,jj,nvc;
604:   PetscScalar       **y;
605:   const PetscScalar *x;

611:   VecGetLocalSize(v,&n);
612:   VecGetLocalSize(s[0],&n2);
613:   VecGetArrayRead(v,&x);
614:   VecGetBlockSize(v,&bs);
615:   if (bs <= 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Input vector does not have a valid blocksize set");

617:   PetscMalloc2(bs,&y,bs,&bss);
618:   nv   = 0;
619:   nvc  = 0;
620:   for (i=0; i<bs; i++) {
621:     VecGetBlockSize(s[i],&bss[i]);
622:     if (bss[i] < 1) bss[i] = 1; /* if user never set it then assume 1  Re: [PETSC #8241] VecStrideGatherAll */
623:     VecGetArray(s[i],&y[i]);
624:     nvc += bss[i];
625:     nv++;
626:     if (nvc > bs)  SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Number of subvectors in subvectors > number of vectors in main vector");
627:     if (nvc == bs) break;
628:   }

630:   n =  n/bs;

632:   jj = 0;
633:   if (addv == INSERT_VALUES) {
634:     for (j=0; j<nv; j++) {
635:       for (k=0; k<bss[j]; k++) {
636:         for (i=0; i<n; i++) y[j][i*bss[j] + k] = x[bs*i+jj+k];
637:       }
638:       jj += bss[j];
639:     }
640:   } else if (addv == ADD_VALUES) {
641:     for (j=0; j<nv; j++) {
642:       for (k=0; k<bss[j]; k++) {
643:         for (i=0; i<n; i++) y[j][i*bss[j] + k] += x[bs*i+jj+k];
644:       }
645:       jj += bss[j];
646:     }
647: #if !defined(PETSC_USE_COMPLEX)
648:   } else if (addv == MAX_VALUES) {
649:     for (j=0; j<nv; j++) {
650:       for (k=0; k<bss[j]; k++) {
651:         for (i=0; i<n; i++) y[j][i*bss[j] + k] = PetscMax(y[j][i*bss[j] + k],x[bs*i+jj+k]);
652:       }
653:       jj += bss[j];
654:     }
655: #endif
656:   } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown insert type");

658:   VecRestoreArrayRead(v,&x);
659:   for (i=0; i<nv; i++) {
660:     VecRestoreArray(s[i],&y[i]);
661:   }

663:   PetscFree2(y,bss);
664:   return(0);
665: }

667: /*@
668:    VecStrideScatterAll - Scatters all the single components from separate vectors into
669:      a multi-component vector.

671:    Collective on Vec

673:    Input Parameters:
674: +  s - the location where the subvectors are stored
675: -  addv - one of ADD_VALUES,INSERT_VALUES,MAX_VALUES

677:    Output Parameter:
678: .  v - the multicomponent vector

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

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

687:    Not optimized; could be easily

689:    Level: advanced

691: .seealso: VecStrideNorm(), VecStrideScatter(), VecStrideMin(), VecStrideMax(), VecStrideGather(),
692:           VecStrideScatterAll()
693: @*/
694: PetscErrorCode  VecStrideScatterAll(Vec s[],Vec v,InsertMode addv)
695: {
696:   PetscErrorCode    ierr;
697:   PetscInt          i,n,n2,bs,j,jj,k,*bss = NULL,nv,nvc;
698:   PetscScalar       *x;
699:   PetscScalar const **y;

705:   VecGetLocalSize(v,&n);
706:   VecGetLocalSize(s[0],&n2);
707:   VecGetArray(v,&x);
708:   VecGetBlockSize(v,&bs);
709:   if (bs <= 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Input vector does not have a valid blocksize set");

711:   PetscMalloc2(bs,(PetscScalar***)&y,bs,&bss);
712:   nv   = 0;
713:   nvc  = 0;
714:   for (i=0; i<bs; i++) {
715:     VecGetBlockSize(s[i],&bss[i]);
716:     if (bss[i] < 1) bss[i] = 1; /* if user never set it then assume 1  Re: [PETSC #8241] VecStrideGatherAll */
717:     VecGetArrayRead(s[i],&y[i]);
718:     nvc += bss[i];
719:     nv++;
720:     if (nvc > bs) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Number of subvectors in subvectors > number of vectors in main vector");
721:     if (nvc == bs) break;
722:   }

724:   n =  n/bs;

726:   jj = 0;
727:   if (addv == INSERT_VALUES) {
728:     for (j=0; j<nv; j++) {
729:       for (k=0; k<bss[j]; k++) {
730:         for (i=0; i<n; i++) x[bs*i+jj+k] = y[j][i*bss[j] + k];
731:       }
732:       jj += bss[j];
733:     }
734:   } else if (addv == ADD_VALUES) {
735:     for (j=0; j<nv; j++) {
736:       for (k=0; k<bss[j]; k++) {
737:         for (i=0; i<n; i++) x[bs*i+jj+k] += y[j][i*bss[j] + k];
738:       }
739:       jj += bss[j];
740:     }
741: #if !defined(PETSC_USE_COMPLEX)
742:   } else if (addv == MAX_VALUES) {
743:     for (j=0; j<nv; j++) {
744:       for (k=0; k<bss[j]; k++) {
745:         for (i=0; i<n; i++) x[bs*i+jj+k] = PetscMax(x[bs*i+jj+k],y[j][i*bss[j] + k]);
746:       }
747:       jj += bss[j];
748:     }
749: #endif
750:   } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown insert type");

752:   VecRestoreArray(v,&x);
753:   for (i=0; i<nv; i++) {
754:     VecRestoreArrayRead(s[i],&y[i]);
755:   }
756:   PetscFree2(*(PetscScalar***)&y,bss);
757:   return(0);
758: }

760: /*@
761:    VecStrideGather - Gathers a single component from a multi-component vector into
762:    another vector.

764:    Collective on Vec

766:    Input Parameters:
767: +  v - the vector
768: .  start - starting point of the subvector (defined by a stride)
769: -  addv - one of ADD_VALUES,INSERT_VALUES,MAX_VALUES

771:    Output Parameter:
772: .  s - the location where the subvector is stored

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

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

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

784:    Not optimized; could be easily

786:    Level: advanced

788: .seealso: VecStrideNorm(), VecStrideScatter(), VecStrideMin(), VecStrideMax(), VecStrideGatherAll(),
789:           VecStrideScatterAll()
790: @*/
791: PetscErrorCode  VecStrideGather(Vec v,PetscInt start,Vec s,InsertMode addv)
792: {

798:   if (start < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative start %D",start);
799:   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);
800:   if (!v->ops->stridegather) SETERRQ(PetscObjectComm((PetscObject)s),PETSC_ERR_SUP,"Not implemented for this Vec class");
801:   (*v->ops->stridegather)(v,start,s,addv);
802:   return(0);
803: }

805: /*@
806:    VecStrideScatter - Scatters a single component from a vector into a multi-component vector.

808:    Collective on Vec

810:    Input Parameters:
811: +  s - the single-component vector
812: .  start - starting point of the subvector (defined by a stride)
813: -  addv - one of ADD_VALUES,INSERT_VALUES,MAX_VALUES

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

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

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

825:    Not optimized; could be easily

827:    Level: advanced

829: .seealso: VecStrideNorm(), VecStrideGather(), VecStrideMin(), VecStrideMax(), VecStrideGatherAll(),
830:           VecStrideScatterAll(), VecStrideSubSetScatter(), VecStrideSubSetGather()
831: @*/
832: PetscErrorCode  VecStrideScatter(Vec s,PetscInt start,Vec v,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->stridescatter) SETERRQ(PetscObjectComm((PetscObject)s),PETSC_ERR_SUP,"Not implemented for this Vec class");
842:   (*v->ops->stridescatter)(s,start,v,addv);
843:   return(0);
844: }

846: /*@
847:    VecStrideSubSetGather - Gathers a subset of components from a multi-component vector into
848:    another vector.

850:    Collective on Vec

852:    Input Parameters:
853: +  v - the vector
854: .  nidx - the number of indices
855: .  idxv - the indices of the components 0 <= idxv[0] ...idxv[nidx-1] < bs(v), they need not be sorted
856: .  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
857: -  addv - one of ADD_VALUES,INSERT_VALUES,MAX_VALUES

859:    Output Parameter:
860: .  s - the location where the subvector is stored

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

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

868:    Not optimized; could be easily

870:    Level: advanced

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

882:   if (nidx == PETSC_DETERMINE) nidx = s->map->bs;
883:   if (!v->ops->stridesubsetgather) SETERRQ(PetscObjectComm((PetscObject)s),PETSC_ERR_SUP,"Not implemented for this Vec class");
884:   (*v->ops->stridesubsetgather)(v,nidx,idxv,idxs,s,addv);
885:   return(0);
886: }

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

891:    Collective on Vec

893:    Input Parameters:
894: +  s - the smaller-component vector
895: .  nidx - the number of indices in idx
896: .  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
897: .  idxv - the indices of the components in the larger-component vector, 0 <= idx[0] ...idx[nidx-1] < bs(v) they need not be sorted
898: -  addv - one of ADD_VALUES,INSERT_VALUES,MAX_VALUES

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

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

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

909:    Not optimized; could be easily

911:    Level: advanced

913: .seealso: VecStrideNorm(), VecStrideGather(), VecStrideGather(), VecStrideSubSetGather(), VecStrideMin(), VecStrideMax(), VecStrideGatherAll(),
914:           VecStrideScatterAll()
915: @*/
916: PetscErrorCode  VecStrideSubSetScatter(Vec s,PetscInt nidx,const PetscInt idxs[],const PetscInt idxv[],Vec v,InsertMode addv)
917: {

923:   if (nidx == PETSC_DETERMINE) nidx = s->map->bs;
924:   if (!v->ops->stridesubsetscatter) SETERRQ(PetscObjectComm((PetscObject)s),PETSC_ERR_SUP,"Not implemented for this Vec class");
925:   (*v->ops->stridesubsetscatter)(s,nidx,idxs,idxv,v,addv);
926:   return(0);
927: }

929: PetscErrorCode  VecStrideGather_Default(Vec v,PetscInt start,Vec s,InsertMode addv)
930: {
932:   PetscInt       i,n,bs,ns;
933:   const PetscScalar *x;
934:   PetscScalar       *y;

937:   VecGetLocalSize(v,&n);
938:   VecGetLocalSize(s,&ns);
939:   VecGetArrayRead(v,&x);
940:   VecGetArray(s,&y);

942:   bs = v->map->bs;
943:   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);
944:   x += start;
945:   n  =  n/bs;

947:   if (addv == INSERT_VALUES) {
948:     for (i=0; i<n; i++) y[i] = x[bs*i];
949:   } else if (addv == ADD_VALUES) {
950:     for (i=0; i<n; i++) y[i] += x[bs*i];
951: #if !defined(PETSC_USE_COMPLEX)
952:   } else if (addv == MAX_VALUES) {
953:     for (i=0; i<n; i++) y[i] = PetscMax(y[i],x[bs*i]);
954: #endif
955:   } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown insert type");

957:   VecRestoreArrayRead(v,&x);
958:   VecRestoreArray(s,&y);
959:   return(0);
960: }

962: PetscErrorCode  VecStrideScatter_Default(Vec s,PetscInt start,Vec v,InsertMode addv)
963: {
964:   PetscErrorCode    ierr;
965:   PetscInt          i,n,bs,ns;
966:   PetscScalar       *x;
967:   const PetscScalar *y;

970:   VecGetLocalSize(v,&n);
971:   VecGetLocalSize(s,&ns);
972:   VecGetArray(v,&x);
973:   VecGetArrayRead(s,&y);

975:   bs = v->map->bs;
976:   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);
977:   x += start;
978:   n  =  n/bs;

980:   if (addv == INSERT_VALUES) {
981:     for (i=0; i<n; i++) x[bs*i] = y[i];
982:   } else if (addv == ADD_VALUES) {
983:     for (i=0; i<n; i++) x[bs*i] += y[i];
984: #if !defined(PETSC_USE_COMPLEX)
985:   } else if (addv == MAX_VALUES) {
986:     for (i=0; i<n; i++) x[bs*i] = PetscMax(y[i],x[bs*i]);
987: #endif
988:   } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown insert type");

990:   VecRestoreArray(v,&x);
991:   VecRestoreArrayRead(s,&y);
992:   return(0);
993: }

995: PetscErrorCode  VecStrideSubSetGather_Default(Vec v,PetscInt nidx,const PetscInt idxv[],const PetscInt idxs[],Vec s,InsertMode addv)
996: {
997:   PetscErrorCode    ierr;
998:   PetscInt          i,j,n,bs,bss,ns;
999:   const PetscScalar *x;
1000:   PetscScalar       *y;

1003:   VecGetLocalSize(v,&n);
1004:   VecGetLocalSize(s,&ns);
1005:   VecGetArrayRead(v,&x);
1006:   VecGetArray(s,&y);

1008:   bs  = v->map->bs;
1009:   bss = s->map->bs;
1010:   n  =  n/bs;

1012:   if (PetscDefined(USE_DEBUG)) {
1013:     if (n != ns/bss) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Incompatible layout of vectors");
1014:     for (j=0; j<nidx; j++) {
1015:       if (idxv[j] < 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"idx[%D] %D is negative",j,idxv[j]);
1016:       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);
1017:     }
1018:     if (!idxs && bss != nidx) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Must provide idxs when not gathering into all locations");
1019:   }

1021:   if (addv == INSERT_VALUES) {
1022:     if (!idxs) {
1023:       for (i=0; i<n; i++) {
1024:         for (j=0; j<bss; j++) y[bss*i+j] = x[bs*i+idxv[j]];
1025:       }
1026:     } else {
1027:       for (i=0; i<n; i++) {
1028:         for (j=0; j<bss; j++) y[bss*i+idxs[j]] = x[bs*i+idxv[j]];
1029:       }
1030:     }
1031:   } else if (addv == ADD_VALUES) {
1032:     if (!idxs) {
1033:       for (i=0; i<n; i++) {
1034:         for (j=0; j<bss; j++) y[bss*i+j] += x[bs*i+idxv[j]];
1035:       }
1036:     } else {
1037:       for (i=0; i<n; i++) {
1038:         for (j=0; j<bss; j++) y[bss*i+idxs[j]] += x[bs*i+idxv[j]];
1039:       }
1040:     }
1041: #if !defined(PETSC_USE_COMPLEX)
1042:   } else if (addv == MAX_VALUES) {
1043:     if (!idxs) {
1044:       for (i=0; i<n; i++) {
1045:         for (j=0; j<bss; j++) y[bss*i+j] = PetscMax(y[bss*i+j],x[bs*i+idxv[j]]);
1046:       }
1047:     } else {
1048:       for (i=0; i<n; i++) {
1049:         for (j=0; j<bss; j++) y[bss*i+idxs[j]] = PetscMax(y[bss*i+idxs[j]],x[bs*i+idxv[j]]);
1050:       }
1051:     }
1052: #endif
1053:   } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown insert type");

1055:   VecRestoreArrayRead(v,&x);
1056:   VecRestoreArray(s,&y);
1057:   return(0);
1058: }

1060: PetscErrorCode  VecStrideSubSetScatter_Default(Vec s,PetscInt nidx,const PetscInt idxs[],const PetscInt idxv[],Vec v,InsertMode addv)
1061: {
1062:   PetscErrorCode    ierr;
1063:   PetscInt          j,i,n,bs,ns,bss;
1064:   PetscScalar       *x;
1065:   const PetscScalar *y;

1068:   VecGetLocalSize(v,&n);
1069:   VecGetLocalSize(s,&ns);
1070:   VecGetArray(v,&x);
1071:   VecGetArrayRead(s,&y);

1073:   bs  = v->map->bs;
1074:   bss = s->map->bs;
1075:   n  =  n/bs;

1077:   if (PetscDefined(USE_DEBUG)) {
1078:     if (n != ns/bss) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Incompatible layout of vectors");
1079:     for (j=0; j<bss; j++) {
1080:       if (idxs) {
1081:         if (idxs[j] < 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"idx[%D] %D is negative",j,idxs[j]);
1082:         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);
1083:       }
1084:     }
1085:     if (!idxs && bss != nidx) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Must provide idxs when not scattering from all locations");
1086:   }

1088:   if (addv == INSERT_VALUES) {
1089:     if (!idxs) {
1090:       for (i=0; i<n; i++) {
1091:         for (j=0; j<bss; j++) x[bs*i + idxv[j]] = y[bss*i+j];
1092:       }
1093:     } else {
1094:       for (i=0; i<n; i++) {
1095:         for (j=0; j<bss; j++) x[bs*i + idxv[j]] = y[bss*i+idxs[j]];
1096:       }
1097:     }
1098:   } else if (addv == ADD_VALUES) {
1099:     if (!idxs) {
1100:       for (i=0; i<n; i++) {
1101:         for (j=0; j<bss; j++) x[bs*i + idxv[j]] += y[bss*i+j];
1102:       }
1103:     } else {
1104:       for (i=0; i<n; i++) {
1105:         for (j=0; j<bss; j++) x[bs*i + idxv[j]] += y[bss*i+idxs[j]];
1106:       }
1107:     }
1108: #if !defined(PETSC_USE_COMPLEX)
1109:   } else if (addv == MAX_VALUES) {
1110:     if (!idxs) {
1111:       for (i=0; i<n; i++) {
1112:         for (j=0; j<bss; j++) x[bs*i + idxv[j]] = PetscMax(y[bss*i+j],x[bs*i + idxv[j]]);
1113:       }
1114:     } else {
1115:       for (i=0; i<n; i++) {
1116:         for (j=0; j<bss; j++) x[bs*i + idxv[j]] = PetscMax(y[bss*i+idxs[j]],x[bs*i + idxv[j]]);
1117:       }
1118:     }
1119: #endif
1120:   } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown insert type");

1122:   VecRestoreArray(v,&x);
1123:   VecRestoreArrayRead(s,&y);
1124:   return(0);
1125: }

1127: PetscErrorCode VecReciprocal_Default(Vec v)
1128: {
1130:   PetscInt       i,n;
1131:   PetscScalar    *x;

1134:   VecGetLocalSize(v,&n);
1135:   VecGetArray(v,&x);
1136:   for (i=0; i<n; i++) {
1137:     if (x[i] != (PetscScalar)0.0) x[i] = (PetscScalar)1.0/x[i];
1138:   }
1139:   VecRestoreArray(v,&x);
1140:   return(0);
1141: }

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

1146:   Not collective

1148:   Input Parameter:
1149: . v - The vector

1151:   Output Parameter:
1152: . v - The vector of exponents

1154:   Level: beginner

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

1158: @*/
1159: PetscErrorCode  VecExp(Vec v)
1160: {
1161:   PetscScalar    *x;
1162:   PetscInt       i, n;

1167:   if (v->ops->exp) {
1168:     (*v->ops->exp)(v);
1169:   } else {
1170:     VecGetLocalSize(v, &n);
1171:     VecGetArray(v, &x);
1172:     for (i = 0; i < n; i++) x[i] = PetscExpScalar(x[i]);
1173:     VecRestoreArray(v, &x);
1174:   }
1175:   return(0);
1176: }

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

1181:   Not collective

1183:   Input Parameter:
1184: . v - The vector

1186:   Output Parameter:
1187: . v - The vector of logs

1189:   Level: beginner

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

1193: @*/
1194: PetscErrorCode  VecLog(Vec v)
1195: {
1196:   PetscScalar    *x;
1197:   PetscInt       i, n;

1202:   if (v->ops->log) {
1203:     (*v->ops->log)(v);
1204:   } else {
1205:     VecGetLocalSize(v, &n);
1206:     VecGetArray(v, &x);
1207:     for (i = 0; i < n; i++) x[i] = PetscLogScalar(x[i]);
1208:     VecRestoreArray(v, &x);
1209:   }
1210:   return(0);
1211: }

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

1216:   Not collective

1218:   Input Parameter:
1219: . v - The vector

1221:   Output Parameter:
1222: . v - The vector square root

1224:   Level: beginner

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

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

1230: @*/
1231: PetscErrorCode  VecSqrtAbs(Vec v)
1232: {
1233:   PetscScalar    *x;
1234:   PetscInt       i, n;

1239:   if (v->ops->sqrt) {
1240:     (*v->ops->sqrt)(v);
1241:   } else {
1242:     VecGetLocalSize(v, &n);
1243:     VecGetArray(v, &x);
1244:     for (i = 0; i < n; i++) x[i] = PetscSqrtReal(PetscAbsScalar(x[i]));
1245:     VecRestoreArray(v, &x);
1246:   }
1247:   return(0);
1248: }

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

1253:   Collective on Vec

1255:   Input Parameters:
1256: + s - first vector
1257: - t - second vector

1259:   Output Parameters:
1260: + dp - s'conj(t)
1261: - nm - t'conj(t)

1263:   Level: advanced

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

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

1270: @*/
1271: PetscErrorCode  VecDotNorm2(Vec s,Vec t,PetscScalar *dp, PetscReal *nm)
1272: {
1273:   const PetscScalar *sx, *tx;
1274:   PetscScalar       dpx = 0.0, nmx = 0.0,work[2],sum[2];
1275:   PetscInt          i, n;
1276:   PetscErrorCode    ierr;

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

1289:   PetscLogEventBegin(VEC_DotNorm2,s,t,0,0);
1290:   if (s->ops->dotnorm2) {
1291:     (*s->ops->dotnorm2)(s,t,dp,&dpx);
1292:     *nm  = PetscRealPart(dpx);
1293:   } else {
1294:     VecGetLocalSize(s, &n);
1295:     VecGetArrayRead(s, &sx);
1296:     VecGetArrayRead(t, &tx);

1298:     for (i = 0; i<n; i++) {
1299:       dpx += sx[i]*PetscConj(tx[i]);
1300:       nmx += tx[i]*PetscConj(tx[i]);
1301:     }
1302:     work[0] = dpx;
1303:     work[1] = nmx;

1305:     MPIU_Allreduce(work,sum,2,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)s));
1306:     *dp  = sum[0];
1307:     *nm  = PetscRealPart(sum[1]);

1309:     VecRestoreArrayRead(t, &tx);
1310:     VecRestoreArrayRead(s, &sx);
1311:     PetscLogFlops(4.0*n);
1312:   }
1313:   PetscLogEventEnd(VEC_DotNorm2,s,t,0,0);
1314:   return(0);
1315: }

1317: /*@
1318:    VecSum - Computes the sum of all the components of a vector.

1320:    Collective on Vec

1322:    Input Parameter:
1323: .  v - the vector

1325:    Output Parameter:
1326: .  sum - the result

1328:    Level: beginner

1330: .seealso: VecMean(), VecNorm()
1331: @*/
1332: PetscErrorCode  VecSum(Vec v,PetscScalar *sum)
1333: {
1334:   PetscErrorCode    ierr;
1335:   PetscInt          i,n;
1336:   const PetscScalar *x;

1341:   *sum = 0.0;
1342:   if (v->ops->sum) {
1343:     (*v->ops->sum)(v,sum);
1344:   } else {
1345:     VecGetLocalSize(v,&n);
1346:     VecGetArrayRead(v,&x);
1347:     for (i=0; i<n; i++) *sum += x[i];
1348:     VecRestoreArrayRead(v,&x);
1349:   }
1350:   MPIU_Allreduce(MPI_IN_PLACE,sum,1,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)v));
1351:   return(0);
1352: }

1354: /*@
1355:    VecMean - Computes the arithmetic mean of all the components of a vector.

1357:    Collective on Vec

1359:    Input Parameter:
1360: .  v - the vector

1362:    Output Parameter:
1363: .  mean - the result

1365:    Level: beginner

1367: .seealso: VecSum(), VecNorm()
1368: @*/
1369: PetscErrorCode  VecMean(Vec v,PetscScalar *mean)
1370: {
1371:   PetscErrorCode    ierr;
1372:   PetscInt          n;

1377:   VecGetSize(v,&n);
1378:   VecSum(v,mean);
1379:   *mean /= n;
1380:   return(0);
1381: }

1383: /*@
1384:    VecImaginaryPart - Replaces a complex vector with its imginary part

1386:    Collective on Vec

1388:    Input Parameter:
1389: .  v - the vector

1391:    Level: beginner

1393: .seealso: VecNorm(), VecRealPart()
1394: @*/
1395: PetscErrorCode  VecImaginaryPart(Vec v)
1396: {
1397:   PetscErrorCode    ierr;
1398:   PetscInt          i,n;
1399:   PetscScalar       *x;

1403:   VecGetLocalSize(v,&n);
1404:   VecGetArray(v,&x);
1405:   for (i=0; i<n; i++) x[i] = PetscImaginaryPart(x[i]);
1406:   VecRestoreArray(v,&x);
1407:   return(0);
1408: }

1410: /*@
1411:    VecRealPart - Replaces a complex vector with its real part

1413:    Collective on Vec

1415:    Input Parameter:
1416: .  v - the vector

1418:    Level: beginner

1420: .seealso: VecNorm(), VecImaginaryPart()
1421: @*/
1422: PetscErrorCode  VecRealPart(Vec v)
1423: {
1424:   PetscErrorCode    ierr;
1425:   PetscInt          i,n;
1426:   PetscScalar       *x;

1430:   VecGetLocalSize(v,&n);
1431:   VecGetArray(v,&x);
1432:   for (i=0; i<n; i++) x[i] = PetscRealPart(x[i]);
1433:   VecRestoreArray(v,&x);
1434:   return(0);
1435: }

1437: /*@
1438:    VecShift - Shifts all of the components of a vector by computing
1439:    x[i] = x[i] + shift.

1441:    Logically Collective on Vec

1443:    Input Parameters:
1444: +  v - the vector
1445: -  shift - the shift

1447:    Level: intermediate

1449: @*/
1450: PetscErrorCode  VecShift(Vec v,PetscScalar shift)
1451: {
1453:   PetscInt       i,n;
1454:   PetscScalar    *x;

1459:   VecSetErrorIfLocked(v,1);
1460:   if (shift == 0.0) return(0);

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

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

1476:    Logically Collective on Vec

1478:    Input Parameters:
1479: .  v - the vector

1481:    Level: intermediate

1483: @*/
1484: PetscErrorCode  VecAbs(Vec v)
1485: {
1487:   PetscInt       i,n;
1488:   PetscScalar    *x;

1492:   VecSetErrorIfLocked(v,1);

1494:   if (v->ops->abs) {
1495:     (*v->ops->abs)(v);
1496:   } else {
1497:     VecGetLocalSize(v,&n);
1498:     VecGetArray(v,&x);
1499:     for (i=0; i<n; i++) x[i] = PetscAbsScalar(x[i]);
1500:     VecRestoreArray(v,&x);
1501:   }
1502:   return(0);
1503: }

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

1508:   Input Parameters:
1509: + vec   - The vector
1510: . order - The ordering
1511: - inv   - The flag for inverting the permutation

1513:   Level: beginner

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

1517: .seealso: MatPermute()
1518: @*/
1519: PetscErrorCode  VecPermute(Vec x, IS row, PetscBool inv)
1520: {
1521:   const PetscScalar *array;
1522:   PetscScalar       *newArray;
1523:   const PetscInt    *idx;
1524:   PetscInt          i,rstart,rend;
1525:   PetscErrorCode    ierr;

1530:   VecSetErrorIfLocked(x,1);
1531:   VecGetOwnershipRange(x,&rstart,&rend);
1532:   ISGetIndices(row, &idx);
1533:   VecGetArrayRead(x, &array);
1534:   PetscMalloc1(x->map->n, &newArray);
1535:   if (PetscDefined(USE_DEBUG)) {
1536:     for (i = 0; i < x->map->n; i++) {
1537:       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]);
1538:     }
1539:   }
1540:   if (!inv) {
1541:     for (i = 0; i < x->map->n; i++) newArray[i] = array[idx[i]-rstart];
1542:   } else {
1543:     for (i = 0; i < x->map->n; i++) newArray[idx[i]-rstart] = array[i];
1544:   }
1545:   VecRestoreArrayRead(x, &array);
1546:   ISRestoreIndices(row, &idx);
1547:   VecReplaceArray(x, newArray);
1548:   return(0);
1549: }

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

1556:    Collective on Vec

1558:    Input Parameters:
1559: +  vec1 - the first vector
1560: -  vec2 - the second vector

1562:    Output Parameter:
1563: .  flg - PETSC_TRUE if the vectors are equal; PETSC_FALSE otherwise.

1565:    Level: intermediate
1566: @*/
1567: PetscErrorCode  VecEqual(Vec vec1,Vec vec2,PetscBool  *flg)
1568: {
1569:   const PetscScalar  *v1,*v2;
1570:   PetscErrorCode     ierr;
1571:   PetscInt           n1,n2,N1,N2;
1572:   PetscBool          flg1;

1578:   if (vec1 == vec2) *flg = PETSC_TRUE;
1579:   else {
1580:     VecGetSize(vec1,&N1);
1581:     VecGetSize(vec2,&N2);
1582:     if (N1 != N2) flg1 = PETSC_FALSE;
1583:     else {
1584:       VecGetLocalSize(vec1,&n1);
1585:       VecGetLocalSize(vec2,&n2);
1586:       if (n1 != n2) flg1 = PETSC_FALSE;
1587:       else {
1588:         VecGetArrayRead(vec1,&v1);
1589:         VecGetArrayRead(vec2,&v2);
1590:         PetscArraycmp(v1,v2,n1,&flg1);
1591:         VecRestoreArrayRead(vec1,&v1);
1592:         VecRestoreArrayRead(vec2,&v2);
1593:       }
1594:     }
1595:     /* combine results from all processors */
1596:     MPIU_Allreduce(&flg1,flg,1,MPIU_BOOL,MPI_MIN,PetscObjectComm((PetscObject)vec1));
1597:   }
1598:   return(0);
1599: }

1601: /*@
1602:    VecUniqueEntries - Compute the number of unique entries, and those entries

1604:    Collective on Vec

1606:    Input Parameter:
1607: .  vec - the vector

1609:    Output Parameters:
1610: +  n - The number of unique entries
1611: -  e - The entries

1613:    Level: intermediate

1615: @*/
1616: PetscErrorCode  VecUniqueEntries(Vec vec, PetscInt *n, PetscScalar **e)
1617: {
1618:   const PetscScalar *v;
1619:   PetscScalar       *tmp, *vals;
1620:   PetscMPIInt       *N, *displs, l;
1621:   PetscInt          ng, m, i, j, p;
1622:   PetscMPIInt       size;
1623:   PetscErrorCode    ierr;

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