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: {
 30:   PetscInt       i,n,bs;
 31:   PetscScalar    *x;

 34:   VecGetLocalSize(v,&n);
 35:   VecGetBlockSize(v,&bs);
 38:   VecGetArray(v,&x);
 39:   for (i=start; i<n; i+=bs) x[i] = s;
 40:   VecRestoreArray(v,&x);
 41:   return 0;
 42: }

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

 48:    Logically Collective on Vec

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

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

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

 61:    Level: advanced

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

 71:   VecGetLocalSize(v,&n);
 72:   VecGetBlockSize(v,&bs);
 75:   VecGetArray(v,&x);
 76:   for (i=start; i<n; i+=bs) x[i] *= scale;
 77:   VecRestoreArray(v,&x);
 78:   return 0;
 79: }

 81: /*@
 82:    VecStrideNorm - Computes the norm of subvector of a vector defined
 83:    by a starting point and a stride.

 85:    Collective on Vec

 87:    Input Parameters:
 88: +  v - the vector
 89: .  start - starting point of the subvector (defined by a stride)
 90: -  ntype - type of norm, one of NORM_1, NORM_2, NORM_INFINITY

 92:    Output Parameter:
 93: .  norm - the norm

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

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

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

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

107:    Level: advanced

109: .seealso: VecNorm(), VecStrideGather(), VecStrideScatter(), VecStrideMin(), VecStrideMax()
110: @*/
111: PetscErrorCode  VecStrideNorm(Vec v,PetscInt start,NormType ntype,PetscReal *nrm)
112: {
113:   PetscInt          i,n,bs;
114:   const PetscScalar *x;
115:   PetscReal         tnorm;

120:   VecGetLocalSize(v,&n);
121:   VecGetBlockSize(v,&bs);
124:   VecGetArrayRead(v,&x);
125:   if (ntype == NORM_2) {
126:     PetscScalar sum = 0.0;
127:     for (i=start; i<n; i+=bs) sum += x[i]*(PetscConj(x[i]));
128:     tnorm = PetscRealPart(sum);
129:     MPIU_Allreduce(&tnorm,nrm,1,MPIU_REAL,MPIU_SUM,PetscObjectComm((PetscObject)v));
130:     *nrm  = PetscSqrtReal(*nrm);
131:   } else if (ntype == NORM_1) {
132:     tnorm = 0.0;
133:     for (i=start; i<n; i+=bs) tnorm += PetscAbsScalar(x[i]);
134:     MPIU_Allreduce(&tnorm,nrm,1,MPIU_REAL,MPIU_SUM,PetscObjectComm((PetscObject)v));
135:   } else if (ntype == NORM_INFINITY) {
136:     tnorm = 0.0;
137:     for (i=start; i<n; i+=bs) {
138:       if (PetscAbsScalar(x[i]) > tnorm) tnorm = PetscAbsScalar(x[i]);
139:     }
140:     MPIU_Allreduce(&tnorm,nrm,1,MPIU_REAL,MPIU_MAX,PetscObjectComm((PetscObject)v));
141:   } else SETERRQ(PetscObjectComm((PetscObject)v),PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown norm type");
142:   VecRestoreArrayRead(v,&x);
143:   return 0;
144: }

146: /*@
147:    VecStrideMax - Computes the maximum of subvector of a vector defined
148:    by a starting point and a stride and optionally its location.

150:    Collective on Vec

152:    Input Parameters:
153: +  v - the vector
154: -  start - starting point of the subvector (defined by a stride)

156:    Output Parameters:
157: +  idex - the location where the maximum occurred  (pass NULL if not required)
158: -  nrm - the maximum value in the subvector

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

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

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

171:    Level: advanced

173: .seealso: VecMax(), VecStrideNorm(), VecStrideGather(), VecStrideScatter(), VecStrideMin()
174: @*/
175: PetscErrorCode  VecStrideMax(Vec v,PetscInt start,PetscInt *idex,PetscReal *nrm)
176: {
177:   PetscInt          i,n,bs,id = -1;
178:   const PetscScalar *x;
179:   PetscReal         max = PETSC_MIN_REAL;

183:   VecGetLocalSize(v,&n);
184:   VecGetBlockSize(v,&bs);
187:   VecGetArrayRead(v,&x);
188:   for (i=start; i<n; i+=bs) {
189:     if (PetscRealPart(x[i]) > max) { max = PetscRealPart(x[i]); id = i;}
190:   }
191:   VecRestoreArrayRead(v,&x);
192: #if defined(PETSC_HAVE_MPIUNI)
193:   *nrm = max;
194:   if (idex) *idex = id;
195: #else
196:   if (!idex) {
197:     MPIU_Allreduce(&max,nrm,1,MPIU_REAL,MPIU_MAX,PetscObjectComm((PetscObject)v));
198:   } else {
199:     struct { PetscReal v; PetscInt i; } in,out;
200:     PetscInt rstart;

202:     VecGetOwnershipRange(v,&rstart,NULL);
203:     in.v  = max;
204:     in.i  = rstart+id;
205:     MPIU_Allreduce(&in,&out,1,MPIU_REAL_INT,MPIU_MAXLOC,PetscObjectComm((PetscObject)v));
206:     *nrm  = out.v;
207:     *idex = out.i;
208:   }
209: #endif
210:   return 0;
211: }

213: /*@
214:    VecStrideMin - Computes the minimum of subvector of a vector defined
215:    by a starting point and a stride and optionally its location.

217:    Collective on Vec

219:    Input Parameters:
220: +  v - the vector
221: -  start - starting point of the subvector (defined by a stride)

223:    Output Parameters:
224: +  idex - the location where the minimum occurred. (pass NULL if not required)
225: -  nrm - the minimum value in the subvector

227:    Level: advanced

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

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

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

240: .seealso: VecMin(), VecStrideNorm(), VecStrideGather(), VecStrideScatter(), VecStrideMax()
241: @*/
242: PetscErrorCode  VecStrideMin(Vec v,PetscInt start,PetscInt *idex,PetscReal *nrm)
243: {
244:   PetscInt          i,n,bs,id = -1;
245:   const PetscScalar *x;
246:   PetscReal         min = PETSC_MAX_REAL;

250:   VecGetLocalSize(v,&n);
251:   VecGetBlockSize(v,&bs);
254:   VecGetArrayRead(v,&x);
255:   for (i=start; i<n; i+=bs) {
256:     if (PetscRealPart(x[i]) < min) { min = PetscRealPart(x[i]); id = i;}
257:   }
258:   VecRestoreArrayRead(v,&x);
259: #if defined(PETSC_HAVE_MPIUNI)
260:   *nrm = min;
261:   if (idex) *idex = id;
262: #else
263:   if (!idex) {
264:     MPIU_Allreduce(&min,nrm,1,MPIU_REAL,MPIU_MIN,PetscObjectComm((PetscObject)v));
265:   } else {
266:     struct { PetscReal v; PetscInt i; } in,out;
267:     PetscInt rstart;

269:     VecGetOwnershipRange(v,&rstart,NULL);
270:     in.v  = min;
271:     in.i  = rstart+id;
272:     MPIU_Allreduce(&in,&out,1,MPIU_REAL_INT,MPIU_MINLOC,PetscObjectComm((PetscObject)v));
273:     *nrm  = out.v;
274:     *idex = out.i;
275:   }
276: #endif
277:   return 0;
278: }

280: /*@
281:    VecStrideScaleAll - Scales the subvectors of a vector defined
282:    by a starting point and a stride.

284:    Logically Collective on Vec

286:    Input Parameters:
287: +  v - the vector
288: -  scales - values to multiply each subvector entry by

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

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

296:    Level: advanced

298: .seealso: VecNorm(), VecStrideScale(), VecScale(), VecStrideGather(), VecStrideScatter(), VecStrideMin(), VecStrideMax()
299: @*/
300: PetscErrorCode  VecStrideScaleAll(Vec v,const PetscScalar *scales)
301: {
302:   PetscInt       i,j,n,bs;
303:   PetscScalar    *x;

307:   VecGetLocalSize(v,&n);
308:   VecGetBlockSize(v,&bs);
309:   VecGetArray(v,&x);
310:   /* need to provide optimized code for each bs */
311:   for (i=0; i<n; i+=bs) {
312:     for (j=0; j<bs; j++) x[i+j] *= scales[j];
313:   }
314:   VecRestoreArray(v,&x);
315:   return 0;
316: }

318: /*@
319:    VecStrideNormAll - Computes the norms of subvectors of a vector defined
320:    by a starting point and a stride.

322:    Collective on Vec

324:    Input Parameters:
325: +  v - the vector
326: -  ntype - type of norm, one of NORM_1, NORM_2, NORM_INFINITY

328:    Output Parameter:
329: .  nrm - the norms

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

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

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

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

342:    Level: advanced

344: .seealso: VecNorm(), VecStrideGather(), VecStrideScatter(), VecStrideMin(), VecStrideMax()
345: @*/
346: PetscErrorCode  VecStrideNormAll(Vec v,NormType ntype,PetscReal nrm[])
347: {
348:   PetscInt          i,j,n,bs;
349:   const PetscScalar *x;
350:   PetscReal         tnorm[128];
351:   MPI_Comm          comm;

356:   VecGetLocalSize(v,&n);
357:   VecGetArrayRead(v,&x);
358:   PetscObjectGetComm((PetscObject)v,&comm);

360:   VecGetBlockSize(v,&bs);

363:   if (ntype == NORM_2) {
364:     PetscScalar sum[128];
365:     for (j=0; j<bs; j++) sum[j] = 0.0;
366:     for (i=0; i<n; i+=bs) {
367:       for (j=0; j<bs; j++) sum[j] += x[i+j]*(PetscConj(x[i+j]));
368:     }
369:     for (j=0; j<bs; j++) tnorm[j]  = PetscRealPart(sum[j]);

371:     MPIU_Allreduce(tnorm,nrm,bs,MPIU_REAL,MPIU_SUM,comm);
372:     for (j=0; j<bs; j++) nrm[j] = PetscSqrtReal(nrm[j]);
373:   } else if (ntype == NORM_1) {
374:     for (j=0; j<bs; j++) tnorm[j] = 0.0;

376:     for (i=0; i<n; i+=bs) {
377:       for (j=0; j<bs; j++) tnorm[j] += PetscAbsScalar(x[i+j]);
378:     }

380:     MPIU_Allreduce(tnorm,nrm,bs,MPIU_REAL,MPIU_SUM,comm);
381:   } else if (ntype == NORM_INFINITY) {
382:     PetscReal tmp;
383:     for (j=0; j<bs; j++) tnorm[j] = 0.0;

385:     for (i=0; i<n; i+=bs) {
386:       for (j=0; j<bs; j++) {
387:         if ((tmp = PetscAbsScalar(x[i+j])) > tnorm[j]) tnorm[j] = tmp;
388:         /* check special case of tmp == NaN */
389:         if (tmp != tmp) {tnorm[j] = tmp; break;}
390:       }
391:     }
392:     MPIU_Allreduce(tnorm,nrm,bs,MPIU_REAL,MPIU_MAX,comm);
393:   } else SETERRQ(comm,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown norm type");
394:   VecRestoreArrayRead(v,&x);
395:   return 0;
396: }

398: /*@
399:    VecStrideMaxAll - Computes the maximums of subvectors of a vector defined
400:    by a starting point and a stride and optionally its location.

402:    Collective on Vec

404:    Input Parameter:
405: .  v - the vector

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

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

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

418:    Level: advanced

420: .seealso: VecMax(), VecStrideNorm(), VecStrideGather(), VecStrideScatter(), VecStrideMin()
421: @*/
422: PetscErrorCode  VecStrideMaxAll(Vec v,PetscInt idex[],PetscReal nrm[])
423: {
424:   PetscInt          i,j,n,bs;
425:   const PetscScalar *x;
426:   PetscReal         max[128],tmp;
427:   MPI_Comm          comm;

432:   VecGetLocalSize(v,&n);
433:   VecGetArrayRead(v,&x);
434:   PetscObjectGetComm((PetscObject)v,&comm);

436:   VecGetBlockSize(v,&bs);

439:   if (!n) {
440:     for (j=0; j<bs; j++) max[j] = PETSC_MIN_REAL;
441:   } else {
442:     for (j=0; j<bs; j++) max[j] = PetscRealPart(x[j]);

444:     for (i=bs; i<n; i+=bs) {
445:       for (j=0; j<bs; j++) {
446:         if ((tmp = PetscRealPart(x[i+j])) > max[j]) max[j] = tmp;
447:       }
448:     }
449:   }
450:   MPIU_Allreduce(max,nrm,bs,MPIU_REAL,MPIU_MAX,comm);

452:   VecRestoreArrayRead(v,&x);
453:   return 0;
454: }

456: /*@
457:    VecStrideMinAll - Computes the minimum of subvector of a vector defined
458:    by a starting point and a stride and optionally its location.

460:    Collective on Vec

462:    Input Parameter:
463: .  v - the vector

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

470:    Level: advanced

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

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

478: .seealso: VecMin(), VecStrideNorm(), VecStrideGather(), VecStrideScatter(), VecStrideMax()
479: @*/
480: PetscErrorCode  VecStrideMinAll(Vec v,PetscInt idex[],PetscReal nrm[])
481: {
482:   PetscInt          i,n,bs,j;
483:   const PetscScalar *x;
484:   PetscReal         min[128],tmp;
485:   MPI_Comm          comm;

490:   VecGetLocalSize(v,&n);
491:   VecGetArrayRead(v,&x);
492:   PetscObjectGetComm((PetscObject)v,&comm);

494:   VecGetBlockSize(v,&bs);

497:   if (!n) {
498:     for (j=0; j<bs; j++) min[j] = PETSC_MAX_REAL;
499:   } else {
500:     for (j=0; j<bs; j++) min[j] = PetscRealPart(x[j]);

502:     for (i=bs; i<n; i+=bs) {
503:       for (j=0; j<bs; j++) {
504:         if ((tmp = PetscRealPart(x[i+j])) < min[j]) min[j] = tmp;
505:       }
506:     }
507:   }
508:   MPIU_Allreduce(min,nrm,bs,MPIU_REAL,MPIU_MIN,comm);

510:   VecRestoreArrayRead(v,&x);
511:   return 0;
512: }

514: /*----------------------------------------------------------------------------------------------*/
515: /*@
516:    VecStrideGatherAll - Gathers all the single components from a multi-component vector into
517:    separate vectors.

519:    Collective on Vec

521:    Input Parameters:
522: +  v - the vector
523: -  addv - one of ADD_VALUES,INSERT_VALUES,MAX_VALUES

525:    Output Parameter:
526: .  s - the location where the subvectors are stored

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

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

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

539:    Not optimized; could be easily

541:    Level: advanced

543: .seealso: VecStrideNorm(), VecStrideScatter(), VecStrideMin(), VecStrideMax(), VecStrideGather(),
544:           VecStrideScatterAll()
545: @*/
546: PetscErrorCode  VecStrideGatherAll(Vec v,Vec s[],InsertMode addv)
547: {
548:   PetscInt          i,n,n2,bs,j,k,*bss = NULL,nv,jj,nvc;
549:   PetscScalar       **y;
550:   const PetscScalar *x;

555:   VecGetLocalSize(v,&n);
556:   VecGetLocalSize(s[0],&n2);
557:   VecGetArrayRead(v,&x);
558:   VecGetBlockSize(v,&bs);

561:   PetscMalloc2(bs,&y,bs,&bss);
562:   nv   = 0;
563:   nvc  = 0;
564:   for (i=0; i<bs; i++) {
565:     VecGetBlockSize(s[i],&bss[i]);
566:     if (bss[i] < 1) bss[i] = 1; /* if user never set it then assume 1  Re: [PETSC #8241] VecStrideGatherAll */
567:     VecGetArray(s[i],&y[i]);
568:     nvc += bss[i];
569:     nv++;
571:     if (nvc == bs) break;
572:   }

574:   n =  n/bs;

576:   jj = 0;
577:   if (addv == INSERT_VALUES) {
578:     for (j=0; j<nv; j++) {
579:       for (k=0; k<bss[j]; k++) {
580:         for (i=0; i<n; i++) y[j][i*bss[j] + k] = x[bs*i+jj+k];
581:       }
582:       jj += bss[j];
583:     }
584:   } else if (addv == ADD_VALUES) {
585:     for (j=0; j<nv; j++) {
586:       for (k=0; k<bss[j]; k++) {
587:         for (i=0; i<n; i++) y[j][i*bss[j] + k] += x[bs*i+jj+k];
588:       }
589:       jj += bss[j];
590:     }
591: #if !defined(PETSC_USE_COMPLEX)
592:   } else if (addv == MAX_VALUES) {
593:     for (j=0; j<nv; j++) {
594:       for (k=0; k<bss[j]; k++) {
595:         for (i=0; i<n; i++) y[j][i*bss[j] + k] = PetscMax(y[j][i*bss[j] + k],x[bs*i+jj+k]);
596:       }
597:       jj += bss[j];
598:     }
599: #endif
600:   } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown insert type");

602:   VecRestoreArrayRead(v,&x);
603:   for (i=0; i<nv; i++) {
604:     VecRestoreArray(s[i],&y[i]);
605:   }

607:   PetscFree2(y,bss);
608:   return 0;
609: }

611: /*@
612:    VecStrideScatterAll - Scatters all the single components from separate vectors into
613:      a multi-component vector.

615:    Collective on Vec

617:    Input Parameters:
618: +  s - the location where the subvectors are stored
619: -  addv - one of ADD_VALUES,INSERT_VALUES,MAX_VALUES

621:    Output Parameter:
622: .  v - the multicomponent vector

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

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

631:    Not optimized; could be easily

633:    Level: advanced

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

647:   VecGetLocalSize(v,&n);
648:   VecGetLocalSize(s[0],&n2);
649:   VecGetArray(v,&x);
650:   VecGetBlockSize(v,&bs);

653:   PetscMalloc2(bs,(PetscScalar***)&y,bs,&bss);
654:   nv   = 0;
655:   nvc  = 0;
656:   for (i=0; i<bs; i++) {
657:     VecGetBlockSize(s[i],&bss[i]);
658:     if (bss[i] < 1) bss[i] = 1; /* if user never set it then assume 1  Re: [PETSC #8241] VecStrideGatherAll */
659:     VecGetArrayRead(s[i],&y[i]);
660:     nvc += bss[i];
661:     nv++;
663:     if (nvc == bs) break;
664:   }

666:   n =  n/bs;

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

694:   VecRestoreArray(v,&x);
695:   for (i=0; i<nv; i++) {
696:     VecRestoreArrayRead(s[i],&y[i]);
697:   }
698:   PetscFree2(*(PetscScalar***)&y,bss);
699:   return 0;
700: }

702: /*@
703:    VecStrideGather - Gathers a single component from a multi-component vector into
704:    another vector.

706:    Collective on Vec

708:    Input Parameters:
709: +  v - the vector
710: .  start - starting point of the subvector (defined by a stride)
711: -  addv - one of ADD_VALUES,INSERT_VALUES,MAX_VALUES

713:    Output Parameter:
714: .  s - the location where the subvector is stored

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

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

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

726:    Not optimized; could be easily

728:    Level: advanced

730: .seealso: VecStrideNorm(), VecStrideScatter(), VecStrideMin(), VecStrideMax(), VecStrideGatherAll(),
731:           VecStrideScatterAll()
732: @*/
733: PetscErrorCode  VecStrideGather(Vec v,PetscInt start,Vec s,InsertMode addv)
734: {
740:   (*v->ops->stridegather)(v,start,s,addv);
741:   return 0;
742: }

744: /*@
745:    VecStrideScatter - Scatters a single component from a vector into a multi-component vector.

747:    Collective on Vec

749:    Input Parameters:
750: +  s - the single-component vector
751: .  start - starting point of the subvector (defined by a stride)
752: -  addv - one of ADD_VALUES,INSERT_VALUES,MAX_VALUES

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

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

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

764:    Not optimized; could be easily

766:    Level: advanced

768: .seealso: VecStrideNorm(), VecStrideGather(), VecStrideMin(), VecStrideMax(), VecStrideGatherAll(),
769:           VecStrideScatterAll(), VecStrideSubSetScatter(), VecStrideSubSetGather()
770: @*/
771: PetscErrorCode  VecStrideScatter(Vec s,PetscInt start,Vec v,InsertMode addv)
772: {
778:   (*v->ops->stridescatter)(s,start,v,addv);
779:   return 0;
780: }

782: /*@
783:    VecStrideSubSetGather - Gathers a subset of components from a multi-component vector into
784:    another vector.

786:    Collective on Vec

788:    Input Parameters:
789: +  v - the vector
790: .  nidx - the number of indices
791: .  idxv - the indices of the components 0 <= idxv[0] ...idxv[nidx-1] < bs(v), they need not be sorted
792: .  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
793: -  addv - one of ADD_VALUES,INSERT_VALUES,MAX_VALUES

795:    Output Parameter:
796: .  s - the location where the subvector is stored

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

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

804:    Not optimized; could be easily

806:    Level: advanced

808: .seealso: VecStrideNorm(), VecStrideScatter(), VecStrideGather(), VecStrideSubSetScatter(), VecStrideMin(), VecStrideMax(), VecStrideGatherAll(),
809:           VecStrideScatterAll()
810: @*/
811: PetscErrorCode  VecStrideSubSetGather(Vec v,PetscInt nidx,const PetscInt idxv[],const PetscInt idxs[],Vec s,InsertMode addv)
812: {
815:   if (nidx == PETSC_DETERMINE) nidx = s->map->bs;
817:   (*v->ops->stridesubsetgather)(v,nidx,idxv,idxs,s,addv);
818:   return 0;
819: }

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

824:    Collective on Vec

826:    Input Parameters:
827: +  s - the smaller-component vector
828: .  nidx - the number of indices in idx
829: .  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
830: .  idxv - the indices of the components in the larger-component vector, 0 <= idx[0] ...idx[nidx-1] < bs(v) they need not be sorted
831: -  addv - one of ADD_VALUES,INSERT_VALUES,MAX_VALUES

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

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

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

842:    Not optimized; could be easily

844:    Level: advanced

846: .seealso: VecStrideNorm(), VecStrideGather(), VecStrideGather(), VecStrideSubSetGather(), VecStrideMin(), VecStrideMax(), VecStrideGatherAll(),
847:           VecStrideScatterAll()
848: @*/
849: PetscErrorCode  VecStrideSubSetScatter(Vec s,PetscInt nidx,const PetscInt idxs[],const PetscInt idxv[],Vec v,InsertMode addv)
850: {
853:   if (nidx == PETSC_DETERMINE) nidx = s->map->bs;
855:   (*v->ops->stridesubsetscatter)(s,nidx,idxs,idxv,v,addv);
856:   return 0;
857: }

859: PetscErrorCode  VecStrideGather_Default(Vec v,PetscInt start,Vec s,InsertMode addv)
860: {
861:   PetscInt       i,n,bs,ns;
862:   const PetscScalar *x;
863:   PetscScalar       *y;

865:   VecGetLocalSize(v,&n);
866:   VecGetLocalSize(s,&ns);
867:   VecGetArrayRead(v,&x);
868:   VecGetArray(s,&y);

870:   bs = v->map->bs;
872:   x += start;
873:   n  =  n/bs;

875:   if (addv == INSERT_VALUES) {
876:     for (i=0; i<n; i++) y[i] = x[bs*i];
877:   } else if (addv == ADD_VALUES) {
878:     for (i=0; i<n; i++) y[i] += x[bs*i];
879: #if !defined(PETSC_USE_COMPLEX)
880:   } else if (addv == MAX_VALUES) {
881:     for (i=0; i<n; i++) y[i] = PetscMax(y[i],x[bs*i]);
882: #endif
883:   } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown insert type");

885:   VecRestoreArrayRead(v,&x);
886:   VecRestoreArray(s,&y);
887:   return 0;
888: }

890: PetscErrorCode  VecStrideScatter_Default(Vec s,PetscInt start,Vec v,InsertMode addv)
891: {
892:   PetscInt          i,n,bs,ns;
893:   PetscScalar       *x;
894:   const PetscScalar *y;

896:   VecGetLocalSize(v,&n);
897:   VecGetLocalSize(s,&ns);
898:   VecGetArray(v,&x);
899:   VecGetArrayRead(s,&y);

901:   bs = v->map->bs;
903:   x += start;
904:   n  =  n/bs;

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

916:   VecRestoreArray(v,&x);
917:   VecRestoreArrayRead(s,&y);
918:   return 0;
919: }

921: PetscErrorCode  VecStrideSubSetGather_Default(Vec v,PetscInt nidx,const PetscInt idxv[],const PetscInt idxs[],Vec s,InsertMode addv)
922: {
923:   PetscInt          i,j,n,bs,bss,ns;
924:   const PetscScalar *x;
925:   PetscScalar       *y;

927:   VecGetLocalSize(v,&n);
928:   VecGetLocalSize(s,&ns);
929:   VecGetArrayRead(v,&x);
930:   VecGetArray(s,&y);

932:   bs  = v->map->bs;
933:   bss = s->map->bs;
934:   n  =  n/bs;

936:   if (PetscDefined(USE_DEBUG)) {
938:     for (j=0; j<nidx; j++) {
941:     }
943:   }

945:   if (addv == INSERT_VALUES) {
946:     if (!idxs) {
947:       for (i=0; i<n; i++) {
948:         for (j=0; j<bss; j++) y[bss*i+j] = x[bs*i+idxv[j]];
949:       }
950:     } else {
951:       for (i=0; i<n; i++) {
952:         for (j=0; j<bss; j++) y[bss*i+idxs[j]] = x[bs*i+idxv[j]];
953:       }
954:     }
955:   } else if (addv == ADD_VALUES) {
956:     if (!idxs) {
957:       for (i=0; i<n; i++) {
958:         for (j=0; j<bss; j++) y[bss*i+j] += x[bs*i+idxv[j]];
959:       }
960:     } else {
961:       for (i=0; i<n; i++) {
962:         for (j=0; j<bss; j++) y[bss*i+idxs[j]] += x[bs*i+idxv[j]];
963:       }
964:     }
965: #if !defined(PETSC_USE_COMPLEX)
966:   } else if (addv == MAX_VALUES) {
967:     if (!idxs) {
968:       for (i=0; i<n; i++) {
969:         for (j=0; j<bss; j++) y[bss*i+j] = PetscMax(y[bss*i+j],x[bs*i+idxv[j]]);
970:       }
971:     } else {
972:       for (i=0; i<n; i++) {
973:         for (j=0; j<bss; j++) y[bss*i+idxs[j]] = PetscMax(y[bss*i+idxs[j]],x[bs*i+idxv[j]]);
974:       }
975:     }
976: #endif
977:   } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown insert type");

979:   VecRestoreArrayRead(v,&x);
980:   VecRestoreArray(s,&y);
981:   return 0;
982: }

984: PetscErrorCode  VecStrideSubSetScatter_Default(Vec s,PetscInt nidx,const PetscInt idxs[],const PetscInt idxv[],Vec v,InsertMode addv)
985: {
986:   PetscInt          j,i,n,bs,ns,bss;
987:   PetscScalar       *x;
988:   const PetscScalar *y;

990:   VecGetLocalSize(v,&n);
991:   VecGetLocalSize(s,&ns);
992:   VecGetArray(v,&x);
993:   VecGetArrayRead(s,&y);

995:   bs  = v->map->bs;
996:   bss = s->map->bs;
997:   n  =  n/bs;

999:   if (PetscDefined(USE_DEBUG)) {
1001:     for (j=0; j<bss; j++) {
1002:       if (idxs) {
1005:       }
1006:     }
1008:   }

1010:   if (addv == INSERT_VALUES) {
1011:     if (!idxs) {
1012:       for (i=0; i<n; i++) {
1013:         for (j=0; j<bss; j++) x[bs*i + idxv[j]] = y[bss*i+j];
1014:       }
1015:     } else {
1016:       for (i=0; i<n; i++) {
1017:         for (j=0; j<bss; j++) x[bs*i + idxv[j]] = y[bss*i+idxs[j]];
1018:       }
1019:     }
1020:   } else if (addv == ADD_VALUES) {
1021:     if (!idxs) {
1022:       for (i=0; i<n; i++) {
1023:         for (j=0; j<bss; j++) x[bs*i + idxv[j]] += y[bss*i+j];
1024:       }
1025:     } else {
1026:       for (i=0; i<n; i++) {
1027:         for (j=0; j<bss; j++) x[bs*i + idxv[j]] += y[bss*i+idxs[j]];
1028:       }
1029:     }
1030: #if !defined(PETSC_USE_COMPLEX)
1031:   } else if (addv == MAX_VALUES) {
1032:     if (!idxs) {
1033:       for (i=0; i<n; i++) {
1034:         for (j=0; j<bss; j++) x[bs*i + idxv[j]] = PetscMax(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++) x[bs*i + idxv[j]] = PetscMax(y[bss*i+idxs[j]],x[bs*i + idxv[j]]);
1039:       }
1040:     }
1041: #endif
1042:   } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown insert type");

1044:   VecRestoreArray(v,&x);
1045:   VecRestoreArrayRead(s,&y);
1046:   return 0;
1047: }

1049: PetscErrorCode VecReciprocal_Default(Vec v)
1050: {
1051:   PetscInt       i,n;
1052:   PetscScalar    *x;

1054:   VecGetLocalSize(v,&n);
1055:   VecGetArray(v,&x);
1056:   for (i=0; i<n; i++) {
1057:     if (x[i] != (PetscScalar)0.0) x[i] = (PetscScalar)1.0/x[i];
1058:   }
1059:   VecRestoreArray(v,&x);
1060:   return 0;
1061: }

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

1066:   Not collective

1068:   Input Parameter:
1069: . v - The vector

1071:   Output Parameter:
1072: . v - The vector of exponents

1074:   Level: beginner

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

1078: @*/
1079: PetscErrorCode  VecExp(Vec v)
1080: {
1081:   PetscScalar    *x;
1082:   PetscInt       i, n;

1085:   if (v->ops->exp) {
1086:     (*v->ops->exp)(v);
1087:   } else {
1088:     VecGetLocalSize(v, &n);
1089:     VecGetArray(v, &x);
1090:     for (i = 0; i < n; i++) x[i] = PetscExpScalar(x[i]);
1091:     VecRestoreArray(v, &x);
1092:   }
1093:   return 0;
1094: }

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

1099:   Not collective

1101:   Input Parameter:
1102: . v - The vector

1104:   Output Parameter:
1105: . v - The vector of logs

1107:   Level: beginner

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

1111: @*/
1112: PetscErrorCode  VecLog(Vec v)
1113: {
1114:   PetscScalar    *x;
1115:   PetscInt       i, n;

1118:   if (v->ops->log) {
1119:     (*v->ops->log)(v);
1120:   } else {
1121:     VecGetLocalSize(v, &n);
1122:     VecGetArray(v, &x);
1123:     for (i = 0; i < n; i++) x[i] = PetscLogScalar(x[i]);
1124:     VecRestoreArray(v, &x);
1125:   }
1126:   return 0;
1127: }

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

1132:   Not collective

1134:   Input Parameter:
1135: . v - The vector

1137:   Output Parameter:
1138: . v - The vector square root

1140:   Level: beginner

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

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

1146: @*/
1147: PetscErrorCode  VecSqrtAbs(Vec v)
1148: {
1149:   PetscScalar    *x;
1150:   PetscInt       i, n;

1153:   if (v->ops->sqrt) {
1154:     (*v->ops->sqrt)(v);
1155:   } else {
1156:     VecGetLocalSize(v, &n);
1157:     VecGetArray(v, &x);
1158:     for (i = 0; i < n; i++) x[i] = PetscSqrtReal(PetscAbsScalar(x[i]));
1159:     VecRestoreArray(v, &x);
1160:   }
1161:   return 0;
1162: }

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

1167:   Collective on Vec

1169:   Input Parameters:
1170: + s - first vector
1171: - t - second vector

1173:   Output Parameters:
1174: + dp - s'conj(t)
1175: - nm - t'conj(t)

1177:   Level: advanced

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

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

1184: @*/
1185: PetscErrorCode  VecDotNorm2(Vec s,Vec t,PetscScalar *dp, PetscReal *nm)
1186: {
1187:   const PetscScalar *sx, *tx;
1188:   PetscScalar       dpx = 0.0, nmx = 0.0,work[2],sum[2];
1189:   PetscInt          i, n;


1201:   PetscLogEventBegin(VEC_DotNorm2,s,t,0,0);
1202:   if (s->ops->dotnorm2) {
1203:     (*s->ops->dotnorm2)(s,t,dp,&dpx);
1204:     *nm  = PetscRealPart(dpx);
1205:   } else {
1206:     VecGetLocalSize(s, &n);
1207:     VecGetArrayRead(s, &sx);
1208:     VecGetArrayRead(t, &tx);

1210:     for (i = 0; i<n; i++) {
1211:       dpx += sx[i]*PetscConj(tx[i]);
1212:       nmx += tx[i]*PetscConj(tx[i]);
1213:     }
1214:     work[0] = dpx;
1215:     work[1] = nmx;

1217:     MPIU_Allreduce(work,sum,2,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)s));
1218:     *dp  = sum[0];
1219:     *nm  = PetscRealPart(sum[1]);

1221:     VecRestoreArrayRead(t, &tx);
1222:     VecRestoreArrayRead(s, &sx);
1223:     PetscLogFlops(4.0*n);
1224:   }
1225:   PetscLogEventEnd(VEC_DotNorm2,s,t,0,0);
1226:   return 0;
1227: }

1229: /*@
1230:    VecSum - Computes the sum of all the components of a vector.

1232:    Collective on Vec

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

1237:    Output Parameter:
1238: .  sum - the result

1240:    Level: beginner

1242: .seealso: VecMean(), VecNorm()
1243: @*/
1244: PetscErrorCode  VecSum(Vec v,PetscScalar *sum)
1245: {
1246:   PetscInt          i,n;
1247:   const PetscScalar *x;

1251:   *sum = 0.0;
1252:   if (v->ops->sum) {
1253:     (*v->ops->sum)(v,sum);
1254:   } else {
1255:     VecGetLocalSize(v,&n);
1256:     VecGetArrayRead(v,&x);
1257:     for (i=0; i<n; i++) *sum += x[i];
1258:     VecRestoreArrayRead(v,&x);
1259:   }
1260:   MPIU_Allreduce(MPI_IN_PLACE,sum,1,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)v));
1261:   return 0;
1262: }

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

1267:    Collective on Vec

1269:    Input Parameter:
1270: .  v - the vector

1272:    Output Parameter:
1273: .  mean - the result

1275:    Level: beginner

1277: .seealso: VecSum(), VecNorm()
1278: @*/
1279: PetscErrorCode  VecMean(Vec v,PetscScalar *mean)
1280: {
1281:   PetscInt          n;

1285:   VecGetSize(v,&n);
1286:   VecSum(v,mean);
1287:   *mean /= n;
1288:   return 0;
1289: }

1291: /*@
1292:    VecImaginaryPart - Replaces a complex vector with its imginary part

1294:    Collective on Vec

1296:    Input Parameter:
1297: .  v - the vector

1299:    Level: beginner

1301: .seealso: VecNorm(), VecRealPart()
1302: @*/
1303: PetscErrorCode  VecImaginaryPart(Vec v)
1304: {
1305:   PetscInt          i,n;
1306:   PetscScalar       *x;

1309:   VecGetLocalSize(v,&n);
1310:   VecGetArray(v,&x);
1311:   for (i=0; i<n; i++) x[i] = PetscImaginaryPart(x[i]);
1312:   VecRestoreArray(v,&x);
1313:   return 0;
1314: }

1316: /*@
1317:    VecRealPart - Replaces a complex vector with its real part

1319:    Collective on Vec

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

1324:    Level: beginner

1326: .seealso: VecNorm(), VecImaginaryPart()
1327: @*/
1328: PetscErrorCode  VecRealPart(Vec v)
1329: {
1330:   PetscInt          i,n;
1331:   PetscScalar       *x;

1334:   VecGetLocalSize(v,&n);
1335:   VecGetArray(v,&x);
1336:   for (i=0; i<n; i++) x[i] = PetscRealPart(x[i]);
1337:   VecRestoreArray(v,&x);
1338:   return 0;
1339: }

1341: /*@
1342:    VecShift - Shifts all of the components of a vector by computing
1343:    x[i] = x[i] + shift.

1345:    Logically Collective on Vec

1347:    Input Parameters:
1348: +  v - the vector
1349: -  shift - the shift

1351:    Level: intermediate

1353: @*/
1354: PetscErrorCode  VecShift(Vec v,PetscScalar shift)
1355: {
1356:   PetscInt       i,n;
1357:   PetscScalar    *x;

1361:   VecSetErrorIfLocked(v,1);
1362:   if (shift == 0.0) return 0;

1364:   if (v->ops->shift) {
1365:     (*v->ops->shift)(v,shift);
1366:   } else {
1367:     VecGetLocalSize(v,&n);
1368:     VecGetArray(v,&x);
1369:     for (i=0; i<n; i++) x[i] += shift;
1370:     VecRestoreArray(v,&x);
1371:   }
1372:   return 0;
1373: }

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

1378:    Logically Collective on Vec

1380:    Input Parameters:
1381: .  v - the vector

1383:    Level: intermediate

1385: @*/
1386: PetscErrorCode  VecAbs(Vec v)
1387: {
1388:   PetscInt       i,n;
1389:   PetscScalar    *x;

1392:   VecSetErrorIfLocked(v,1);

1394:   if (v->ops->abs) {
1395:     (*v->ops->abs)(v);
1396:   } else {
1397:     VecGetLocalSize(v,&n);
1398:     VecGetArray(v,&x);
1399:     for (i=0; i<n; i++) x[i] = PetscAbsScalar(x[i]);
1400:     VecRestoreArray(v,&x);
1401:   }
1402:   return 0;
1403: }

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

1408:   Input Parameters:
1409: + vec   - The vector
1410: . order - The ordering
1411: - inv   - The flag for inverting the permutation

1413:   Level: beginner

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

1417: .seealso: MatPermute()
1418: @*/
1419: PetscErrorCode  VecPermute(Vec x, IS row, PetscBool inv)
1420: {
1421:   const PetscScalar *array;
1422:   PetscScalar       *newArray;
1423:   const PetscInt    *idx;
1424:   PetscInt          i,rstart,rend;

1428:   VecSetErrorIfLocked(x,1);
1429:   VecGetOwnershipRange(x,&rstart,&rend);
1430:   ISGetIndices(row, &idx);
1431:   VecGetArrayRead(x, &array);
1432:   PetscMalloc1(x->map->n, &newArray);
1433:   if (PetscDefined(USE_DEBUG)) {
1434:     for (i = 0; i < x->map->n; i++) {
1436:     }
1437:   }
1438:   if (!inv) {
1439:     for (i = 0; i < x->map->n; i++) newArray[i] = array[idx[i]-rstart];
1440:   } else {
1441:     for (i = 0; i < x->map->n; i++) newArray[idx[i]-rstart] = array[i];
1442:   }
1443:   VecRestoreArrayRead(x, &array);
1444:   ISRestoreIndices(row, &idx);
1445:   VecReplaceArray(x, newArray);
1446:   return 0;
1447: }

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

1454:    Collective on Vec

1456:    Input Parameters:
1457: +  vec1 - the first vector
1458: -  vec2 - the second vector

1460:    Output Parameter:
1461: .  flg - PETSC_TRUE if the vectors are equal; PETSC_FALSE otherwise.

1463:    Level: intermediate
1464: @*/
1465: PetscErrorCode  VecEqual(Vec vec1,Vec vec2,PetscBool  *flg)
1466: {
1467:   const PetscScalar  *v1,*v2;
1468:   PetscInt           n1,n2,N1,N2;
1469:   PetscBool          flg1;

1474:   if (vec1 == vec2) *flg = PETSC_TRUE;
1475:   else {
1476:     VecGetSize(vec1,&N1);
1477:     VecGetSize(vec2,&N2);
1478:     if (N1 != N2) flg1 = PETSC_FALSE;
1479:     else {
1480:       VecGetLocalSize(vec1,&n1);
1481:       VecGetLocalSize(vec2,&n2);
1482:       if (n1 != n2) flg1 = PETSC_FALSE;
1483:       else {
1484:         VecGetArrayRead(vec1,&v1);
1485:         VecGetArrayRead(vec2,&v2);
1486:         PetscArraycmp(v1,v2,n1,&flg1);
1487:         VecRestoreArrayRead(vec1,&v1);
1488:         VecRestoreArrayRead(vec2,&v2);
1489:       }
1490:     }
1491:     /* combine results from all processors */
1492:     MPIU_Allreduce(&flg1,flg,1,MPIU_BOOL,MPI_MIN,PetscObjectComm((PetscObject)vec1));
1493:   }
1494:   return 0;
1495: }

1497: /*@
1498:    VecUniqueEntries - Compute the number of unique entries, and those entries

1500:    Collective on Vec

1502:    Input Parameter:
1503: .  vec - the vector

1505:    Output Parameters:
1506: +  n - The number of unique entries
1507: -  e - The entries

1509:    Level: intermediate

1511: @*/
1512: PetscErrorCode  VecUniqueEntries(Vec vec, PetscInt *n, PetscScalar **e)
1513: {
1514:   const PetscScalar *v;
1515:   PetscScalar       *tmp, *vals;
1516:   PetscMPIInt       *N, *displs, l;
1517:   PetscInt          ng, m, i, j, p;
1518:   PetscMPIInt       size;

1522:   MPI_Comm_size(PetscObjectComm((PetscObject) vec), &size);
1523:   VecGetLocalSize(vec, &m);
1524:   VecGetArrayRead(vec, &v);
1525:   PetscMalloc2(m,&tmp,size,&N);
1526:   for (i = 0, j = 0, l = 0; i < m; ++i) {
1527:     /* Can speed this up with sorting */
1528:     for (j = 0; j < l; ++j) {
1529:       if (v[i] == tmp[j]) break;
1530:     }
1531:     if (j == l) {
1532:       tmp[j] = v[i];
1533:       ++l;
1534:     }
1535:   }
1536:   VecRestoreArrayRead(vec, &v);
1537:   /* Gather serial results */
1538:   MPI_Allgather(&l, 1, MPI_INT, N, 1, MPI_INT, PetscObjectComm((PetscObject) vec));
1539:   for (p = 0, ng = 0; p < size; ++p) {
1540:     ng += N[p];
1541:   }
1542:   PetscMalloc2(ng,&vals,size+1,&displs);
1543:   for (p = 1, displs[0] = 0; p <= size; ++p) {
1544:     displs[p] = displs[p-1] + N[p-1];
1545:   }
1546:   MPI_Allgatherv(tmp, l, MPIU_SCALAR, vals, N, displs, MPIU_SCALAR, PetscObjectComm((PetscObject) vec));
1547:   /* Find unique entries */
1548: #ifdef PETSC_USE_COMPLEX
1549:   SETERRQ(PetscObjectComm((PetscObject) vec), PETSC_ERR_SUP, "Does not work with complex numbers");
1550: #else
1551:   *n = displs[size];
1552:   PetscSortRemoveDupsReal(n, (PetscReal *) vals);
1553:   if (e) {
1555:     PetscMalloc1(*n, e);
1556:     for (i = 0; i < *n; ++i) {
1557:       (*e)[i] = vals[i];
1558:     }
1559:   }
1560:   PetscFree2(vals,displs);
1561:   PetscFree2(tmp,N);
1562:   return 0;
1563: #endif
1564: }