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;

 36:   VecGetLocalSize(v,&n);
 37:   VecGetBlockSize(v,&bs);
 38:   if (start < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative start %D",start);
 39:   else if (start >= bs) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Start of stride subvector (%D) is too large for stride %D",start,bs);
 40:   VecGetArray(v,&x);
 41:   for (i=start; i<n; i+=bs) x[i] = s;
 42:   VecRestoreArray(v,&x);
 43:   return(0);
 44: }

 46: /*@
 47:    VecStrideScale - Scales a subvector of a vector defined
 48:    by a starting point and a stride.

 50:    Logically Collective on Vec

 52:    Input Parameters:
 53: +  v - the vector
 54: .  start - starting point of the subvector (defined by a stride)
 55: -  scale - value to multiply each subvector entry by

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

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

 63:    Level: advanced

 65: .seealso: VecNorm(), VecStrideGather(), VecStrideScatter(), VecStrideMin(), VecStrideMax(), VecStrideScale()
 66: @*/
 67: PetscErrorCode  VecStrideScale(Vec v,PetscInt start,PetscScalar scale)
 68: {
 70:   PetscInt       i,n,bs;
 71:   PetscScalar    *x;

 75:   VecGetLocalSize(v,&n);
 76:   VecGetBlockSize(v,&bs);
 77:   if (start < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative start %D",start);
 78:   else if (start >= bs) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Start of stride subvector (%D) is too large for stride %D",start,bs);
 79:   VecGetArray(v,&x);
 80:   for (i=start; i<n; i+=bs) x[i] *= scale;
 81:   VecRestoreArray(v,&x);
 82:   return(0);
 83: }

 85: /*@
 86:    VecStrideNorm - Computes the norm of subvector of a vector defined
 87:    by a starting point and a stride.

 89:    Collective on Vec

 91:    Input Parameters:
 92: +  v - the vector
 93: .  start - starting point of the subvector (defined by a stride)
 94: -  ntype - type of norm, one of NORM_1, NORM_2, NORM_INFINITY

 96:    Output Parameter:
 97: .  norm - the norm

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

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

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

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

111:    Level: advanced

113: .seealso: VecNorm(), VecStrideGather(), VecStrideScatter(), VecStrideMin(), VecStrideMax()
114: @*/
115: PetscErrorCode  VecStrideNorm(Vec v,PetscInt start,NormType ntype,PetscReal *nrm)
116: {
117:   PetscErrorCode    ierr;
118:   PetscInt          i,n,bs;
119:   const PetscScalar *x;
120:   PetscReal         tnorm;

126:   VecGetLocalSize(v,&n);
127:   VecGetBlockSize(v,&bs);
128:   if (start < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative start %D",start);
129:   else if (start >= bs) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Start of stride subvector (%D) is too large for stride %D",start,bs);
130:   VecGetArrayRead(v,&x);
131:   if (ntype == NORM_2) {
132:     PetscScalar sum = 0.0;
133:     for (i=start; i<n; i+=bs) sum += x[i]*(PetscConj(x[i]));
134:     tnorm = PetscRealPart(sum);
135:     MPIU_Allreduce(&tnorm,nrm,1,MPIU_REAL,MPIU_SUM,PetscObjectComm((PetscObject)v));
136:     *nrm  = PetscSqrtReal(*nrm);
137:   } else if (ntype == NORM_1) {
138:     tnorm = 0.0;
139:     for (i=start; i<n; i+=bs) tnorm += PetscAbsScalar(x[i]);
140:     MPIU_Allreduce(&tnorm,nrm,1,MPIU_REAL,MPIU_SUM,PetscObjectComm((PetscObject)v));
141:   } else if (ntype == NORM_INFINITY) {
142:     tnorm = 0.0;
143:     for (i=start; i<n; i+=bs) {
144:       if (PetscAbsScalar(x[i]) > tnorm) tnorm = PetscAbsScalar(x[i]);
145:     }
146:     MPIU_Allreduce(&tnorm,nrm,1,MPIU_REAL,MPIU_MAX,PetscObjectComm((PetscObject)v));
147:   } else SETERRQ(PetscObjectComm((PetscObject)v),PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown norm type");
148:   VecRestoreArrayRead(v,&x);
149:   return(0);
150: }

152: /*@
153:    VecStrideMax - Computes the maximum of subvector of a vector defined
154:    by a starting point and a stride and optionally its location.

156:    Collective on Vec

158:    Input Parameters:
159: +  v - the vector
160: -  start - starting point of the subvector (defined by a stride)

162:    Output Parameters:
163: +  idex - the location where the maximum occurred  (pass NULL if not required)
164: -  nrm - the maximum value in the subvector

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

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

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

177:    Level: advanced

179: .seealso: VecMax(), VecStrideNorm(), VecStrideGather(), VecStrideScatter(), VecStrideMin()
180: @*/
181: PetscErrorCode  VecStrideMax(Vec v,PetscInt start,PetscInt *idex,PetscReal *nrm)
182: {
183:   PetscErrorCode    ierr;
184:   PetscInt          i,n,bs,id = -1;
185:   const PetscScalar *x;
186:   PetscReal         max = PETSC_MIN_REAL;

191:   VecGetLocalSize(v,&n);
192:   VecGetBlockSize(v,&bs);
193:   if (start < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative start %D",start);
194:   else if (start >= bs) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Start of stride subvector (%D) is too large for stride %D",start,bs);
195:   VecGetArrayRead(v,&x);
196:   for (i=start; i<n; i+=bs) {
197:     if (PetscRealPart(x[i]) > max) { max = PetscRealPart(x[i]); id = i;}
198:   }
199:   VecRestoreArrayRead(v,&x);
200: #if defined(PETSC_HAVE_MPIUNI)
201:   *nrm = max;
202:   if (idex) *idex = id;
203: #else
204:   if (!idex) {
205:     MPIU_Allreduce(&max,nrm,1,MPIU_REAL,MPIU_MAX,PetscObjectComm((PetscObject)v));
206:   } else {
207:     struct { PetscReal v; PetscInt i; } in,out;
208:     PetscInt rstart;

210:     VecGetOwnershipRange(v,&rstart,NULL);
211:     in.v  = max;
212:     in.i  = rstart+id;
213:     MPIU_Allreduce(&in,&out,1,MPIU_REAL_INT,MPIU_MAXLOC,PetscObjectComm((PetscObject)v));
214:     *nrm  = out.v;
215:     *idex = out.i;
216:   }
217: #endif
218:   return(0);
219: }

221: /*@
222:    VecStrideMin - Computes the minimum of subvector of a vector defined
223:    by a starting point and a stride and optionally its location.

225:    Collective on Vec

227:    Input Parameters:
228: +  v - the vector
229: -  start - starting point of the subvector (defined by a stride)

231:    Output Parameters:
232: +  idex - the location where the minimum occurred. (pass NULL if not required)
233: -  nrm - the minimum value in the subvector

235:    Level: advanced

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

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

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

248: .seealso: VecMin(), VecStrideNorm(), VecStrideGather(), VecStrideScatter(), VecStrideMax()
249: @*/
250: PetscErrorCode  VecStrideMin(Vec v,PetscInt start,PetscInt *idex,PetscReal *nrm)
251: {
252:   PetscErrorCode    ierr;
253:   PetscInt          i,n,bs,id = -1;
254:   const PetscScalar *x;
255:   PetscReal         min = PETSC_MAX_REAL;

260:   VecGetLocalSize(v,&n);
261:   VecGetBlockSize(v,&bs);
262:   if (start < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative start %D",start);
263:   else if (start >= bs) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Start of stride subvector (%D) is too large for stride %D",start,bs);
264:   VecGetArrayRead(v,&x);
265:   for (i=start; i<n; i+=bs) {
266:     if (PetscRealPart(x[i]) < min) { min = PetscRealPart(x[i]); id = i;}
267:   }
268:   VecRestoreArrayRead(v,&x);
269: #if defined(PETSC_HAVE_MPIUNI)
270:   *nrm = min;
271:   if (idex) *idex = id;
272: #else
273:   if (!idex) {
274:     MPIU_Allreduce(&min,nrm,1,MPIU_REAL,MPIU_MIN,PetscObjectComm((PetscObject)v));
275:   } else {
276:     struct { PetscReal v; PetscInt i; } in,out;
277:     PetscInt rstart;

279:     VecGetOwnershipRange(v,&rstart,NULL);
280:     in.v  = min;
281:     in.i  = rstart+id;
282:     MPIU_Allreduce(&in,&out,1,MPIU_REAL_INT,MPIU_MINLOC,PetscObjectComm((PetscObject)v));
283:     *nrm  = out.v;
284:     *idex = out.i;
285:   }
286: #endif
287:   return(0);
288: }

290: /*@
291:    VecStrideScaleAll - Scales the subvectors of a vector defined
292:    by a starting point and a stride.

294:    Logically Collective on Vec

296:    Input Parameters:
297: +  v - the vector
298: -  scales - values to multiply each subvector entry by

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

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

306:    Level: advanced

308: .seealso: VecNorm(), VecStrideScale(), VecScale(), VecStrideGather(), VecStrideScatter(), VecStrideMin(), VecStrideMax()
309: @*/
310: PetscErrorCode  VecStrideScaleAll(Vec v,const PetscScalar *scales)
311: {
313:   PetscInt       i,j,n,bs;
314:   PetscScalar    *x;

319:   VecGetLocalSize(v,&n);
320:   VecGetBlockSize(v,&bs);
321:   VecGetArray(v,&x);
322:   /* need to provide optimized code for each bs */
323:   for (i=0; i<n; i+=bs) {
324:     for (j=0; j<bs; j++) x[i+j] *= scales[j];
325:   }
326:   VecRestoreArray(v,&x);
327:   return(0);
328: }

330: /*@
331:    VecStrideNormAll - Computes the norms of subvectors of a vector defined
332:    by a starting point and a stride.

334:    Collective on Vec

336:    Input Parameters:
337: +  v - the vector
338: -  ntype - type of norm, one of NORM_1, NORM_2, NORM_INFINITY

340:    Output Parameter:
341: .  nrm - the norms

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:    If x is the array representing the vector x then this computes the norm
348:    of the array (x[start],x[start+stride],x[start+2*stride], ....) for each start < stride

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

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

354:    Level: advanced

356: .seealso: VecNorm(), VecStrideGather(), VecStrideScatter(), VecStrideMin(), VecStrideMax()
357: @*/
358: PetscErrorCode  VecStrideNormAll(Vec v,NormType ntype,PetscReal nrm[])
359: {
360:   PetscErrorCode    ierr;
361:   PetscInt          i,j,n,bs;
362:   const PetscScalar *x;
363:   PetscReal         tnorm[128];
364:   MPI_Comm          comm;

370:   VecGetLocalSize(v,&n);
371:   VecGetArrayRead(v,&x);
372:   PetscObjectGetComm((PetscObject)v,&comm);

374:   VecGetBlockSize(v,&bs);
375:   if (bs > 128) SETERRQ(comm,PETSC_ERR_SUP,"Currently supports only blocksize up to 128");

377:   if (ntype == NORM_2) {
378:     PetscScalar sum[128];
379:     for (j=0; j<bs; j++) sum[j] = 0.0;
380:     for (i=0; i<n; i+=bs) {
381:       for (j=0; j<bs; j++) sum[j] += x[i+j]*(PetscConj(x[i+j]));
382:     }
383:     for (j=0; j<bs; j++) tnorm[j]  = PetscRealPart(sum[j]);

385:     MPIU_Allreduce(tnorm,nrm,bs,MPIU_REAL,MPIU_SUM,comm);
386:     for (j=0; j<bs; j++) nrm[j] = PetscSqrtReal(nrm[j]);
387:   } else if (ntype == NORM_1) {
388:     for (j=0; j<bs; j++) tnorm[j] = 0.0;

390:     for (i=0; i<n; i+=bs) {
391:       for (j=0; j<bs; j++) tnorm[j] += PetscAbsScalar(x[i+j]);
392:     }

394:     MPIU_Allreduce(tnorm,nrm,bs,MPIU_REAL,MPIU_SUM,comm);
395:   } else if (ntype == NORM_INFINITY) {
396:     PetscReal tmp;
397:     for (j=0; j<bs; j++) tnorm[j] = 0.0;

399:     for (i=0; i<n; i+=bs) {
400:       for (j=0; j<bs; j++) {
401:         if ((tmp = PetscAbsScalar(x[i+j])) > tnorm[j]) tnorm[j] = tmp;
402:         /* check special case of tmp == NaN */
403:         if (tmp != tmp) {tnorm[j] = tmp; break;}
404:       }
405:     }
406:     MPIU_Allreduce(tnorm,nrm,bs,MPIU_REAL,MPIU_MAX,comm);
407:   } else SETERRQ(comm,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown norm type");
408:   VecRestoreArrayRead(v,&x);
409:   return(0);
410: }

412: /*@
413:    VecStrideMaxAll - Computes the maximums of subvectors of a vector defined
414:    by a starting point and a stride and optionally its location.

416:    Collective on Vec

418:    Input Parameter:
419: .  v - the vector

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

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

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

432:    Level: advanced

434: .seealso: VecMax(), VecStrideNorm(), VecStrideGather(), VecStrideScatter(), VecStrideMin()
435: @*/
436: PetscErrorCode  VecStrideMaxAll(Vec v,PetscInt idex[],PetscReal nrm[])
437: {
438:   PetscErrorCode    ierr;
439:   PetscInt          i,j,n,bs;
440:   const PetscScalar *x;
441:   PetscReal         max[128],tmp;
442:   MPI_Comm          comm;

447:   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");
448:   VecGetLocalSize(v,&n);
449:   VecGetArrayRead(v,&x);
450:   PetscObjectGetComm((PetscObject)v,&comm);

452:   VecGetBlockSize(v,&bs);
453:   if (bs > 128) SETERRQ(comm,PETSC_ERR_SUP,"Currently supports only blocksize up to 128");

455:   if (!n) {
456:     for (j=0; j<bs; j++) max[j] = PETSC_MIN_REAL;
457:   } else {
458:     for (j=0; j<bs; j++) max[j] = PetscRealPart(x[j]);

460:     for (i=bs; i<n; i+=bs) {
461:       for (j=0; j<bs; j++) {
462:         if ((tmp = PetscRealPart(x[i+j])) > max[j]) max[j] = tmp;
463:       }
464:     }
465:   }
466:   MPIU_Allreduce(max,nrm,bs,MPIU_REAL,MPIU_MAX,comm);

468:   VecRestoreArrayRead(v,&x);
469:   return(0);
470: }

472: /*@
473:    VecStrideMinAll - Computes the minimum of subvector of a vector defined
474:    by a starting point and a stride and optionally its location.

476:    Collective on Vec

478:    Input Parameter:
479: .  v - the vector

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

486:    Level: advanced

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

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

494: .seealso: VecMin(), VecStrideNorm(), VecStrideGather(), VecStrideScatter(), VecStrideMax()
495: @*/
496: PetscErrorCode  VecStrideMinAll(Vec v,PetscInt idex[],PetscReal nrm[])
497: {
498:   PetscErrorCode    ierr;
499:   PetscInt          i,n,bs,j;
500:   const PetscScalar *x;
501:   PetscReal         min[128],tmp;
502:   MPI_Comm          comm;

507:   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");
508:   VecGetLocalSize(v,&n);
509:   VecGetArrayRead(v,&x);
510:   PetscObjectGetComm((PetscObject)v,&comm);

512:   VecGetBlockSize(v,&bs);
513:   if (bs > 128) SETERRQ(comm,PETSC_ERR_SUP,"Currently supports only blocksize up to 128");

515:   if (!n) {
516:     for (j=0; j<bs; j++) min[j] = PETSC_MAX_REAL;
517:   } else {
518:     for (j=0; j<bs; j++) min[j] = PetscRealPart(x[j]);

520:     for (i=bs; i<n; i+=bs) {
521:       for (j=0; j<bs; j++) {
522:         if ((tmp = PetscRealPart(x[i+j])) < min[j]) min[j] = tmp;
523:       }
524:     }
525:   }
526:   MPIU_Allreduce(min,nrm,bs,MPIU_REAL,MPIU_MIN,comm);

528:   VecRestoreArrayRead(v,&x);
529:   return(0);
530: }

532: /*----------------------------------------------------------------------------------------------*/
533: /*@
534:    VecStrideGatherAll - Gathers all the single components from a multi-component vector into
535:    separate vectors.

537:    Collective on Vec

539:    Input Parameters:
540: +  v - the vector
541: -  addv - one of ADD_VALUES,INSERT_VALUES,MAX_VALUES

543:    Output Parameter:
544: .  s - the location where the subvectors are stored

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

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

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

557:    Not optimized; could be easily

559:    Level: advanced

561: .seealso: VecStrideNorm(), VecStrideScatter(), VecStrideMin(), VecStrideMax(), VecStrideGather(),
562:           VecStrideScatterAll()
563: @*/
564: PetscErrorCode  VecStrideGatherAll(Vec v,Vec s[],InsertMode addv)
565: {
566:   PetscErrorCode    ierr;
567:   PetscInt          i,n,n2,bs,j,k,*bss = NULL,nv,jj,nvc;
568:   PetscScalar       **y;
569:   const PetscScalar *x;

575:   VecGetLocalSize(v,&n);
576:   VecGetLocalSize(s[0],&n2);
577:   VecGetArrayRead(v,&x);
578:   VecGetBlockSize(v,&bs);
579:   if (bs <= 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Input vector does not have a valid blocksize set");

581:   PetscMalloc2(bs,&y,bs,&bss);
582:   nv   = 0;
583:   nvc  = 0;
584:   for (i=0; i<bs; i++) {
585:     VecGetBlockSize(s[i],&bss[i]);
586:     if (bss[i] < 1) bss[i] = 1; /* if user never set it then assume 1  Re: [PETSC #8241] VecStrideGatherAll */
587:     VecGetArray(s[i],&y[i]);
588:     nvc += bss[i];
589:     nv++;
590:     if (nvc > bs)  SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Number of subvectors in subvectors > number of vectors in main vector");
591:     if (nvc == bs) break;
592:   }

594:   n =  n/bs;

596:   jj = 0;
597:   if (addv == INSERT_VALUES) {
598:     for (j=0; j<nv; j++) {
599:       for (k=0; k<bss[j]; k++) {
600:         for (i=0; i<n; i++) y[j][i*bss[j] + k] = x[bs*i+jj+k];
601:       }
602:       jj += bss[j];
603:     }
604:   } else if (addv == ADD_VALUES) {
605:     for (j=0; j<nv; j++) {
606:       for (k=0; k<bss[j]; k++) {
607:         for (i=0; i<n; i++) y[j][i*bss[j] + k] += x[bs*i+jj+k];
608:       }
609:       jj += bss[j];
610:     }
611: #if !defined(PETSC_USE_COMPLEX)
612:   } else if (addv == MAX_VALUES) {
613:     for (j=0; j<nv; j++) {
614:       for (k=0; k<bss[j]; k++) {
615:         for (i=0; i<n; i++) y[j][i*bss[j] + k] = PetscMax(y[j][i*bss[j] + k],x[bs*i+jj+k]);
616:       }
617:       jj += bss[j];
618:     }
619: #endif
620:   } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown insert type");

622:   VecRestoreArrayRead(v,&x);
623:   for (i=0; i<nv; i++) {
624:     VecRestoreArray(s[i],&y[i]);
625:   }

627:   PetscFree2(y,bss);
628:   return(0);
629: }

631: /*@
632:    VecStrideScatterAll - Scatters all the single components from separate vectors into
633:      a multi-component vector.

635:    Collective on Vec

637:    Input Parameters:
638: +  s - the location where the subvectors are stored
639: -  addv - one of ADD_VALUES,INSERT_VALUES,MAX_VALUES

641:    Output Parameter:
642: .  v - the multicomponent vector

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

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

651:    Not optimized; could be easily

653:    Level: advanced

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

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

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

688:   n =  n/bs;

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

716:   VecRestoreArray(v,&x);
717:   for (i=0; i<nv; i++) {
718:     VecRestoreArrayRead(s[i],&y[i]);
719:   }
720:   PetscFree2(*(PetscScalar***)&y,bss);
721:   return(0);
722: }

724: /*@
725:    VecStrideGather - Gathers a single component from a multi-component vector into
726:    another vector.

728:    Collective on Vec

730:    Input Parameters:
731: +  v - the vector
732: .  start - starting point of the subvector (defined by a stride)
733: -  addv - one of ADD_VALUES,INSERT_VALUES,MAX_VALUES

735:    Output Parameter:
736: .  s - the location where the subvector is stored

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

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

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

748:    Not optimized; could be easily

750:    Level: advanced

752: .seealso: VecStrideNorm(), VecStrideScatter(), VecStrideMin(), VecStrideMax(), VecStrideGatherAll(),
753:           VecStrideScatterAll()
754: @*/
755: PetscErrorCode  VecStrideGather(Vec v,PetscInt start,Vec s,InsertMode addv)
756: {

762:   if (start < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative start %D",start);
763:   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);
764:   if (!v->ops->stridegather) SETERRQ(PetscObjectComm((PetscObject)s),PETSC_ERR_SUP,"Not implemented for this Vec class");
765:   (*v->ops->stridegather)(v,start,s,addv);
766:   return(0);
767: }

769: /*@
770:    VecStrideScatter - Scatters a single component from a vector into a multi-component vector.

772:    Collective on Vec

774:    Input Parameters:
775: +  s - the single-component vector
776: .  start - starting point of the subvector (defined by a stride)
777: -  addv - one of ADD_VALUES,INSERT_VALUES,MAX_VALUES

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

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

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

789:    Not optimized; could be easily

791:    Level: advanced

793: .seealso: VecStrideNorm(), VecStrideGather(), VecStrideMin(), VecStrideMax(), VecStrideGatherAll(),
794:           VecStrideScatterAll(), VecStrideSubSetScatter(), VecStrideSubSetGather()
795: @*/
796: PetscErrorCode  VecStrideScatter(Vec s,PetscInt start,Vec v,InsertMode addv)
797: {

803:   if (start < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative start %D",start);
804:   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);
805:   if (!v->ops->stridescatter) SETERRQ(PetscObjectComm((PetscObject)s),PETSC_ERR_SUP,"Not implemented for this Vec class");
806:   (*v->ops->stridescatter)(s,start,v,addv);
807:   return(0);
808: }

810: /*@
811:    VecStrideSubSetGather - Gathers a subset of components from a multi-component vector into
812:    another vector.

814:    Collective on Vec

816:    Input Parameters:
817: +  v - the vector
818: .  nidx - the number of indices
819: .  idxv - the indices of the components 0 <= idxv[0] ...idxv[nidx-1] < bs(v), they need not be sorted
820: .  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
821: -  addv - one of ADD_VALUES,INSERT_VALUES,MAX_VALUES

823:    Output Parameter:
824: .  s - the location where the subvector is stored

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

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

832:    Not optimized; could be easily

834:    Level: advanced

836: .seealso: VecStrideNorm(), VecStrideScatter(), VecStrideGather(), VecStrideSubSetScatter(), VecStrideMin(), VecStrideMax(), VecStrideGatherAll(),
837:           VecStrideScatterAll()
838: @*/
839: PetscErrorCode  VecStrideSubSetGather(Vec v,PetscInt nidx,const PetscInt idxv[],const PetscInt idxs[],Vec s,InsertMode addv)
840: {

846:   if (nidx == PETSC_DETERMINE) nidx = s->map->bs;
847:   if (!v->ops->stridesubsetgather) SETERRQ(PetscObjectComm((PetscObject)s),PETSC_ERR_SUP,"Not implemented for this Vec class");
848:   (*v->ops->stridesubsetgather)(v,nidx,idxv,idxs,s,addv);
849:   return(0);
850: }

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

855:    Collective on Vec

857:    Input Parameters:
858: +  s - the smaller-component vector
859: .  nidx - the number of indices in idx
860: .  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
861: .  idxv - the indices of the components in the larger-component vector, 0 <= idx[0] ...idx[nidx-1] < bs(v) they need not be sorted
862: -  addv - one of ADD_VALUES,INSERT_VALUES,MAX_VALUES

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

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

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

873:    Not optimized; could be easily

875:    Level: advanced

877: .seealso: VecStrideNorm(), VecStrideGather(), VecStrideGather(), VecStrideSubSetGather(), VecStrideMin(), VecStrideMax(), VecStrideGatherAll(),
878:           VecStrideScatterAll()
879: @*/
880: PetscErrorCode  VecStrideSubSetScatter(Vec s,PetscInt nidx,const PetscInt idxs[],const PetscInt idxv[],Vec v,InsertMode addv)
881: {

887:   if (nidx == PETSC_DETERMINE) nidx = s->map->bs;
888:   if (!v->ops->stridesubsetscatter) SETERRQ(PetscObjectComm((PetscObject)s),PETSC_ERR_SUP,"Not implemented for this Vec class");
889:   (*v->ops->stridesubsetscatter)(s,nidx,idxs,idxv,v,addv);
890:   return(0);
891: }

893: PetscErrorCode  VecStrideGather_Default(Vec v,PetscInt start,Vec s,InsertMode addv)
894: {
896:   PetscInt       i,n,bs,ns;
897:   const PetscScalar *x;
898:   PetscScalar       *y;

901:   VecGetLocalSize(v,&n);
902:   VecGetLocalSize(s,&ns);
903:   VecGetArrayRead(v,&x);
904:   VecGetArray(s,&y);

906:   bs = v->map->bs;
907:   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);
908:   x += start;
909:   n  =  n/bs;

911:   if (addv == INSERT_VALUES) {
912:     for (i=0; i<n; i++) y[i] = x[bs*i];
913:   } else if (addv == ADD_VALUES) {
914:     for (i=0; i<n; i++) y[i] += x[bs*i];
915: #if !defined(PETSC_USE_COMPLEX)
916:   } else if (addv == MAX_VALUES) {
917:     for (i=0; i<n; i++) y[i] = PetscMax(y[i],x[bs*i]);
918: #endif
919:   } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown insert type");

921:   VecRestoreArrayRead(v,&x);
922:   VecRestoreArray(s,&y);
923:   return(0);
924: }

926: PetscErrorCode  VecStrideScatter_Default(Vec s,PetscInt start,Vec v,InsertMode addv)
927: {
928:   PetscErrorCode    ierr;
929:   PetscInt          i,n,bs,ns;
930:   PetscScalar       *x;
931:   const PetscScalar *y;

934:   VecGetLocalSize(v,&n);
935:   VecGetLocalSize(s,&ns);
936:   VecGetArray(v,&x);
937:   VecGetArrayRead(s,&y);

939:   bs = v->map->bs;
940:   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);
941:   x += start;
942:   n  =  n/bs;

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

954:   VecRestoreArray(v,&x);
955:   VecRestoreArrayRead(s,&y);
956:   return(0);
957: }

959: PetscErrorCode  VecStrideSubSetGather_Default(Vec v,PetscInt nidx,const PetscInt idxv[],const PetscInt idxs[],Vec s,InsertMode addv)
960: {
961:   PetscErrorCode    ierr;
962:   PetscInt          i,j,n,bs,bss,ns;
963:   const PetscScalar *x;
964:   PetscScalar       *y;

967:   VecGetLocalSize(v,&n);
968:   VecGetLocalSize(s,&ns);
969:   VecGetArrayRead(v,&x);
970:   VecGetArray(s,&y);

972:   bs  = v->map->bs;
973:   bss = s->map->bs;
974:   n  =  n/bs;

976:   if (PetscDefined(USE_DEBUG)) {
977:     if (n != ns/bss) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Incompatible layout of vectors");
978:     for (j=0; j<nidx; j++) {
979:       if (idxv[j] < 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"idx[%D] %D is negative",j,idxv[j]);
980:       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);
981:     }
982:     if (!idxs && bss != nidx) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Must provide idxs when not gathering into all locations");
983:   }

985:   if (addv == INSERT_VALUES) {
986:     if (!idxs) {
987:       for (i=0; i<n; i++) {
988:         for (j=0; j<bss; j++) y[bss*i+j] = x[bs*i+idxv[j]];
989:       }
990:     } else {
991:       for (i=0; i<n; i++) {
992:         for (j=0; j<bss; j++) y[bss*i+idxs[j]] = x[bs*i+idxv[j]];
993:       }
994:     }
995:   } else if (addv == ADD_VALUES) {
996:     if (!idxs) {
997:       for (i=0; i<n; i++) {
998:         for (j=0; j<bss; j++) y[bss*i+j] += x[bs*i+idxv[j]];
999:       }
1000:     } else {
1001:       for (i=0; i<n; i++) {
1002:         for (j=0; j<bss; j++) y[bss*i+idxs[j]] += x[bs*i+idxv[j]];
1003:       }
1004:     }
1005: #if !defined(PETSC_USE_COMPLEX)
1006:   } else if (addv == MAX_VALUES) {
1007:     if (!idxs) {
1008:       for (i=0; i<n; i++) {
1009:         for (j=0; j<bss; j++) y[bss*i+j] = PetscMax(y[bss*i+j],x[bs*i+idxv[j]]);
1010:       }
1011:     } else {
1012:       for (i=0; i<n; i++) {
1013:         for (j=0; j<bss; j++) y[bss*i+idxs[j]] = PetscMax(y[bss*i+idxs[j]],x[bs*i+idxv[j]]);
1014:       }
1015:     }
1016: #endif
1017:   } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown insert type");

1019:   VecRestoreArrayRead(v,&x);
1020:   VecRestoreArray(s,&y);
1021:   return(0);
1022: }

1024: PetscErrorCode  VecStrideSubSetScatter_Default(Vec s,PetscInt nidx,const PetscInt idxs[],const PetscInt idxv[],Vec v,InsertMode addv)
1025: {
1026:   PetscErrorCode    ierr;
1027:   PetscInt          j,i,n,bs,ns,bss;
1028:   PetscScalar       *x;
1029:   const PetscScalar *y;

1032:   VecGetLocalSize(v,&n);
1033:   VecGetLocalSize(s,&ns);
1034:   VecGetArray(v,&x);
1035:   VecGetArrayRead(s,&y);

1037:   bs  = v->map->bs;
1038:   bss = s->map->bs;
1039:   n  =  n/bs;

1041:   if (PetscDefined(USE_DEBUG)) {
1042:     if (n != ns/bss) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Incompatible layout of vectors");
1043:     for (j=0; j<bss; j++) {
1044:       if (idxs) {
1045:         if (idxs[j] < 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"idx[%D] %D is negative",j,idxs[j]);
1046:         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);
1047:       }
1048:     }
1049:     if (!idxs && bss != nidx) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Must provide idxs when not scattering from all locations");
1050:   }

1052:   if (addv == INSERT_VALUES) {
1053:     if (!idxs) {
1054:       for (i=0; i<n; i++) {
1055:         for (j=0; j<bss; j++) x[bs*i + idxv[j]] = y[bss*i+j];
1056:       }
1057:     } else {
1058:       for (i=0; i<n; i++) {
1059:         for (j=0; j<bss; j++) x[bs*i + idxv[j]] = y[bss*i+idxs[j]];
1060:       }
1061:     }
1062:   } else if (addv == ADD_VALUES) {
1063:     if (!idxs) {
1064:       for (i=0; i<n; i++) {
1065:         for (j=0; j<bss; j++) x[bs*i + idxv[j]] += y[bss*i+j];
1066:       }
1067:     } else {
1068:       for (i=0; i<n; i++) {
1069:         for (j=0; j<bss; j++) x[bs*i + idxv[j]] += y[bss*i+idxs[j]];
1070:       }
1071:     }
1072: #if !defined(PETSC_USE_COMPLEX)
1073:   } else if (addv == MAX_VALUES) {
1074:     if (!idxs) {
1075:       for (i=0; i<n; i++) {
1076:         for (j=0; j<bss; j++) x[bs*i + idxv[j]] = PetscMax(y[bss*i+j],x[bs*i + idxv[j]]);
1077:       }
1078:     } else {
1079:       for (i=0; i<n; i++) {
1080:         for (j=0; j<bss; j++) x[bs*i + idxv[j]] = PetscMax(y[bss*i+idxs[j]],x[bs*i + idxv[j]]);
1081:       }
1082:     }
1083: #endif
1084:   } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown insert type");

1086:   VecRestoreArray(v,&x);
1087:   VecRestoreArrayRead(s,&y);
1088:   return(0);
1089: }

1091: PetscErrorCode VecReciprocal_Default(Vec v)
1092: {
1094:   PetscInt       i,n;
1095:   PetscScalar    *x;

1098:   VecGetLocalSize(v,&n);
1099:   VecGetArray(v,&x);
1100:   for (i=0; i<n; i++) {
1101:     if (x[i] != (PetscScalar)0.0) x[i] = (PetscScalar)1.0/x[i];
1102:   }
1103:   VecRestoreArray(v,&x);
1104:   return(0);
1105: }

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

1110:   Not collective

1112:   Input Parameter:
1113: . v - The vector

1115:   Output Parameter:
1116: . v - The vector of exponents

1118:   Level: beginner

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

1122: @*/
1123: PetscErrorCode  VecExp(Vec v)
1124: {
1125:   PetscScalar    *x;
1126:   PetscInt       i, n;

1131:   if (v->ops->exp) {
1132:     (*v->ops->exp)(v);
1133:   } else {
1134:     VecGetLocalSize(v, &n);
1135:     VecGetArray(v, &x);
1136:     for (i = 0; i < n; i++) x[i] = PetscExpScalar(x[i]);
1137:     VecRestoreArray(v, &x);
1138:   }
1139:   return(0);
1140: }

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

1145:   Not collective

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

1150:   Output Parameter:
1151: . v - The vector of logs

1153:   Level: beginner

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

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

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

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

1180:   Not collective

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

1185:   Output Parameter:
1186: . v - The vector square root

1188:   Level: beginner

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

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

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

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

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

1217:   Collective on Vec

1219:   Input Parameters:
1220: + s - first vector
1221: - t - second vector

1223:   Output Parameters:
1224: + dp - s'conj(t)
1225: - nm - t'conj(t)

1227:   Level: advanced

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

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

1234: @*/
1235: PetscErrorCode  VecDotNorm2(Vec s,Vec t,PetscScalar *dp, PetscReal *nm)
1236: {
1237:   const PetscScalar *sx, *tx;
1238:   PetscScalar       dpx = 0.0, nmx = 0.0,work[2],sum[2];
1239:   PetscInt          i, n;
1240:   PetscErrorCode    ierr;

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

1253:   PetscLogEventBegin(VEC_DotNorm2,s,t,0,0);
1254:   if (s->ops->dotnorm2) {
1255:     (*s->ops->dotnorm2)(s,t,dp,&dpx);
1256:     *nm  = PetscRealPart(dpx);
1257:   } else {
1258:     VecGetLocalSize(s, &n);
1259:     VecGetArrayRead(s, &sx);
1260:     VecGetArrayRead(t, &tx);

1262:     for (i = 0; i<n; i++) {
1263:       dpx += sx[i]*PetscConj(tx[i]);
1264:       nmx += tx[i]*PetscConj(tx[i]);
1265:     }
1266:     work[0] = dpx;
1267:     work[1] = nmx;

1269:     MPIU_Allreduce(work,sum,2,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)s));
1270:     *dp  = sum[0];
1271:     *nm  = PetscRealPart(sum[1]);

1273:     VecRestoreArrayRead(t, &tx);
1274:     VecRestoreArrayRead(s, &sx);
1275:     PetscLogFlops(4.0*n);
1276:   }
1277:   PetscLogEventEnd(VEC_DotNorm2,s,t,0,0);
1278:   return(0);
1279: }

1281: /*@
1282:    VecSum - Computes the sum of all the components of a vector.

1284:    Collective on Vec

1286:    Input Parameter:
1287: .  v - the vector

1289:    Output Parameter:
1290: .  sum - the result

1292:    Level: beginner

1294: .seealso: VecMean(), VecNorm()
1295: @*/
1296: PetscErrorCode  VecSum(Vec v,PetscScalar *sum)
1297: {
1298:   PetscErrorCode    ierr;
1299:   PetscInt          i,n;
1300:   const PetscScalar *x;

1305:   *sum = 0.0;
1306:   if (v->ops->sum) {
1307:     (*v->ops->sum)(v,sum);
1308:   } else {
1309:     VecGetLocalSize(v,&n);
1310:     VecGetArrayRead(v,&x);
1311:     for (i=0; i<n; i++) *sum += x[i];
1312:     VecRestoreArrayRead(v,&x);
1313:   }
1314:   MPIU_Allreduce(MPI_IN_PLACE,sum,1,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)v));
1315:   return(0);
1316: }

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

1321:    Collective on Vec

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

1326:    Output Parameter:
1327: .  mean - the result

1329:    Level: beginner

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

1341:   VecGetSize(v,&n);
1342:   VecSum(v,mean);
1343:   *mean /= n;
1344:   return(0);
1345: }

1347: /*@
1348:    VecImaginaryPart - Replaces a complex vector with its imginary part

1350:    Collective on Vec

1352:    Input Parameter:
1353: .  v - the vector

1355:    Level: beginner

1357: .seealso: VecNorm(), VecRealPart()
1358: @*/
1359: PetscErrorCode  VecImaginaryPart(Vec v)
1360: {
1361:   PetscErrorCode    ierr;
1362:   PetscInt          i,n;
1363:   PetscScalar       *x;

1367:   VecGetLocalSize(v,&n);
1368:   VecGetArray(v,&x);
1369:   for (i=0; i<n; i++) x[i] = PetscImaginaryPart(x[i]);
1370:   VecRestoreArray(v,&x);
1371:   return(0);
1372: }

1374: /*@
1375:    VecRealPart - Replaces a complex vector with its real part

1377:    Collective on Vec

1379:    Input Parameter:
1380: .  v - the vector

1382:    Level: beginner

1384: .seealso: VecNorm(), VecImaginaryPart()
1385: @*/
1386: PetscErrorCode  VecRealPart(Vec v)
1387: {
1388:   PetscErrorCode    ierr;
1389:   PetscInt          i,n;
1390:   PetscScalar       *x;

1394:   VecGetLocalSize(v,&n);
1395:   VecGetArray(v,&x);
1396:   for (i=0; i<n; i++) x[i] = PetscRealPart(x[i]);
1397:   VecRestoreArray(v,&x);
1398:   return(0);
1399: }

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

1405:    Logically Collective on Vec

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

1411:    Level: intermediate

1413: @*/
1414: PetscErrorCode  VecShift(Vec v,PetscScalar shift)
1415: {
1417:   PetscInt       i,n;
1418:   PetscScalar    *x;

1423:   VecSetErrorIfLocked(v,1);
1424:   if (shift == 0.0) return(0);

1426:   if (v->ops->shift) {
1427:     (*v->ops->shift)(v,shift);
1428:   } else {
1429:     VecGetLocalSize(v,&n);
1430:     VecGetArray(v,&x);
1431:     for (i=0; i<n; i++) x[i] += shift;
1432:     VecRestoreArray(v,&x);
1433:   }
1434:   return(0);
1435: }

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

1440:    Logically Collective on Vec

1442:    Input Parameters:
1443: .  v - the vector

1445:    Level: intermediate

1447: @*/
1448: PetscErrorCode  VecAbs(Vec v)
1449: {
1451:   PetscInt       i,n;
1452:   PetscScalar    *x;

1456:   VecSetErrorIfLocked(v,1);

1458:   if (v->ops->abs) {
1459:     (*v->ops->abs)(v);
1460:   } else {
1461:     VecGetLocalSize(v,&n);
1462:     VecGetArray(v,&x);
1463:     for (i=0; i<n; i++) x[i] = PetscAbsScalar(x[i]);
1464:     VecRestoreArray(v,&x);
1465:   }
1466:   return(0);
1467: }

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

1472:   Input Parameters:
1473: + vec   - The vector
1474: . order - The ordering
1475: - inv   - The flag for inverting the permutation

1477:   Level: beginner

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

1481: .seealso: MatPermute()
1482: @*/
1483: PetscErrorCode  VecPermute(Vec x, IS row, PetscBool inv)
1484: {
1485:   const PetscScalar *array;
1486:   PetscScalar       *newArray;
1487:   const PetscInt    *idx;
1488:   PetscInt          i,rstart,rend;
1489:   PetscErrorCode    ierr;

1494:   VecSetErrorIfLocked(x,1);
1495:   VecGetOwnershipRange(x,&rstart,&rend);
1496:   ISGetIndices(row, &idx);
1497:   VecGetArrayRead(x, &array);
1498:   PetscMalloc1(x->map->n, &newArray);
1499:   if (PetscDefined(USE_DEBUG)) {
1500:     for (i = 0; i < x->map->n; i++) {
1501:       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]);
1502:     }
1503:   }
1504:   if (!inv) {
1505:     for (i = 0; i < x->map->n; i++) newArray[i] = array[idx[i]-rstart];
1506:   } else {
1507:     for (i = 0; i < x->map->n; i++) newArray[idx[i]-rstart] = array[i];
1508:   }
1509:   VecRestoreArrayRead(x, &array);
1510:   ISRestoreIndices(row, &idx);
1511:   VecReplaceArray(x, newArray);
1512:   return(0);
1513: }

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

1520:    Collective on Vec

1522:    Input Parameters:
1523: +  vec1 - the first vector
1524: -  vec2 - the second vector

1526:    Output Parameter:
1527: .  flg - PETSC_TRUE if the vectors are equal; PETSC_FALSE otherwise.

1529:    Level: intermediate
1530: @*/
1531: PetscErrorCode  VecEqual(Vec vec1,Vec vec2,PetscBool  *flg)
1532: {
1533:   const PetscScalar  *v1,*v2;
1534:   PetscErrorCode     ierr;
1535:   PetscInt           n1,n2,N1,N2;
1536:   PetscBool          flg1;

1542:   if (vec1 == vec2) *flg = PETSC_TRUE;
1543:   else {
1544:     VecGetSize(vec1,&N1);
1545:     VecGetSize(vec2,&N2);
1546:     if (N1 != N2) flg1 = PETSC_FALSE;
1547:     else {
1548:       VecGetLocalSize(vec1,&n1);
1549:       VecGetLocalSize(vec2,&n2);
1550:       if (n1 != n2) flg1 = PETSC_FALSE;
1551:       else {
1552:         VecGetArrayRead(vec1,&v1);
1553:         VecGetArrayRead(vec2,&v2);
1554:         PetscArraycmp(v1,v2,n1,&flg1);
1555:         VecRestoreArrayRead(vec1,&v1);
1556:         VecRestoreArrayRead(vec2,&v2);
1557:       }
1558:     }
1559:     /* combine results from all processors */
1560:     MPIU_Allreduce(&flg1,flg,1,MPIU_BOOL,MPI_MIN,PetscObjectComm((PetscObject)vec1));
1561:   }
1562:   return(0);
1563: }

1565: /*@
1566:    VecUniqueEntries - Compute the number of unique entries, and those entries

1568:    Collective on Vec

1570:    Input Parameter:
1571: .  vec - the vector

1573:    Output Parameters:
1574: +  n - The number of unique entries
1575: -  e - The entries

1577:    Level: intermediate

1579: @*/
1580: PetscErrorCode  VecUniqueEntries(Vec vec, PetscInt *n, PetscScalar **e)
1581: {
1582:   const PetscScalar *v;
1583:   PetscScalar       *tmp, *vals;
1584:   PetscMPIInt       *N, *displs, l;
1585:   PetscInt          ng, m, i, j, p;
1586:   PetscMPIInt       size;
1587:   PetscErrorCode    ierr;

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