Actual source code: pbvec.c

petsc-3.4.4 2014-03-13
  2: /*
  3:    This file contains routines for Parallel vector operations.
  4:  */
  5: #include <../src/vec/vec/impls/mpi/pvecimpl.h>   /*I  "petscvec.h"   I*/

  9: static PetscErrorCode VecDot_MPI(Vec xin,Vec yin,PetscScalar *z)
 10: {
 11:   PetscScalar    sum,work;

 15:   VecDot_Seq(xin,yin,&work);
 16:   MPI_Allreduce(&work,&sum,1,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)xin));
 17:   *z   = sum;
 18:   return(0);
 19: }

 23: static PetscErrorCode VecTDot_MPI(Vec xin,Vec yin,PetscScalar *z)
 24: {
 25:   PetscScalar    sum,work;

 29:   VecTDot_Seq(xin,yin,&work);
 30:   MPI_Allreduce(&work,&sum,1,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)xin));
 31:   *z   = sum;
 32:   return(0);
 33: }

 35: extern PetscErrorCode VecView_MPI_Draw(Vec,PetscViewer);

 39: static PetscErrorCode VecPlaceArray_MPI(Vec vin,const PetscScalar *a)
 40: {
 42:   Vec_MPI        *v = (Vec_MPI*)vin->data;

 45:   if (v->unplacedarray) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"VecPlaceArray() was already called on this vector, without a call to VecResetArray()");
 46:   v->unplacedarray = v->array;  /* save previous array so reset can bring it back */
 47:   v->array         = (PetscScalar*)a;
 48:   if (v->localrep) {
 49:     VecPlaceArray(v->localrep,a);
 50:   }
 51:   return(0);
 52: }

 54: extern PetscErrorCode VecGetValues_MPI(Vec,PetscInt,const PetscInt[],PetscScalar[]);

 58: static PetscErrorCode VecDuplicate_MPI(Vec win,Vec *v)
 59: {
 61:   Vec_MPI        *vw,*w = (Vec_MPI*)win->data;
 62:   PetscScalar    *array;

 65:   VecCreate(PetscObjectComm((PetscObject)win),v);
 66:   PetscLayoutReference(win->map,&(*v)->map);

 68:   VecCreate_MPI_Private(*v,PETSC_TRUE,w->nghost,0);
 69:   vw   = (Vec_MPI*)(*v)->data;
 70:   PetscMemcpy((*v)->ops,win->ops,sizeof(struct _VecOps));

 72:   /* save local representation of the parallel vector (and scatter) if it exists */
 73:   if (w->localrep) {
 74:     VecGetArray(*v,&array);
 75:     VecCreateSeqWithArray(PETSC_COMM_SELF,win->map->bs,win->map->n+w->nghost,array,&vw->localrep);
 76:     PetscMemcpy(vw->localrep->ops,w->localrep->ops,sizeof(struct _VecOps));
 77:     VecRestoreArray(*v,&array);
 78:     PetscLogObjectParent(*v,vw->localrep);

 80:     vw->localupdate = w->localupdate;
 81:     if (vw->localupdate) {
 82:       PetscObjectReference((PetscObject)vw->localupdate);
 83:     }
 84:   }

 86:   /* New vector should inherit stashing property of parent */
 87:   (*v)->stash.donotstash   = win->stash.donotstash;
 88:   (*v)->stash.ignorenegidx = win->stash.ignorenegidx;

 90:   PetscObjectListDuplicate(((PetscObject)win)->olist,&((PetscObject)(*v))->olist);
 91:   PetscFunctionListDuplicate(((PetscObject)win)->qlist,&((PetscObject)(*v))->qlist);

 93:   (*v)->map->bs   = win->map->bs;
 94:   (*v)->bstash.bs = win->bstash.bs;
 95:   return(0);
 96: }

 98: extern PetscErrorCode VecSetOption_MPI(Vec,VecOption,PetscBool);
 99: extern PetscErrorCode VecResetArray_MPI(Vec);

101: static struct _VecOps DvOps = { VecDuplicate_MPI, /* 1 */
102:                                 VecDuplicateVecs_Default,
103:                                 VecDestroyVecs_Default,
104:                                 VecDot_MPI,
105:                                 VecMDot_MPI,
106:                                 VecNorm_MPI,
107:                                 VecTDot_MPI,
108:                                 VecMTDot_MPI,
109:                                 VecScale_Seq,
110:                                 VecCopy_Seq, /* 10 */
111:                                 VecSet_Seq,
112:                                 VecSwap_Seq,
113:                                 VecAXPY_Seq,
114:                                 VecAXPBY_Seq,
115:                                 VecMAXPY_Seq,
116:                                 VecAYPX_Seq,
117:                                 VecWAXPY_Seq,
118:                                 VecAXPBYPCZ_Seq,
119:                                 VecPointwiseMult_Seq,
120:                                 VecPointwiseDivide_Seq,
121:                                 VecSetValues_MPI, /* 20 */
122:                                 VecAssemblyBegin_MPI,
123:                                 VecAssemblyEnd_MPI,
124:                                 0,
125:                                 VecGetSize_MPI,
126:                                 VecGetSize_Seq,
127:                                 0,
128:                                 VecMax_MPI,
129:                                 VecMin_MPI,
130:                                 VecSetRandom_Seq,
131:                                 VecSetOption_MPI,
132:                                 VecSetValuesBlocked_MPI,
133:                                 VecDestroy_MPI,
134:                                 VecView_MPI,
135:                                 VecPlaceArray_MPI,
136:                                 VecReplaceArray_Seq,
137:                                 VecDot_Seq,
138:                                 VecTDot_Seq,
139:                                 VecNorm_Seq,
140:                                 VecMDot_Seq,
141:                                 VecMTDot_Seq,
142:                                 VecLoad_Default,
143:                                 VecReciprocal_Default,
144:                                 VecConjugate_Seq,
145:                                 0,
146:                                 0,
147:                                 VecResetArray_MPI,
148:                                 0,
149:                                 VecMaxPointwiseDivide_Seq,
150:                                 VecPointwiseMax_Seq,
151:                                 VecPointwiseMaxAbs_Seq,
152:                                 VecPointwiseMin_Seq,
153:                                 VecGetValues_MPI,
154:                                 0,
155:                                 0,
156:                                 0,
157:                                 0,
158:                                 0,
159:                                 0,
160:                                 VecStrideGather_Default,
161:                                 VecStrideScatter_Default,
162:                                 0,
163:                                 0,
164:                                 0,
165:                                 0,
166:                                 0
167: };

171: /*
172:     VecCreate_MPI_Private - Basic create routine called by VecCreate_MPI() (i.e. VecCreateMPI()),
173:     VecCreateMPIWithArray(), VecCreate_Shared() (i.e. VecCreateShared()), VecCreateGhost(),
174:     VecDuplicate_MPI(), VecCreateGhostWithArray(), VecDuplicate_MPI(), and VecDuplicate_Shared()

176:     If alloc is true and array is NULL then this routine allocates the space, otherwise
177:     no space is allocated.
178: */
179: PetscErrorCode VecCreate_MPI_Private(Vec v,PetscBool alloc,PetscInt nghost,const PetscScalar array[])
180: {
181:   Vec_MPI        *s;

185:   PetscNewLog(v,Vec_MPI,&s);
186:   v->data        = (void*)s;
187:   PetscMemcpy(v->ops,&DvOps,sizeof(DvOps));
188:   s->nghost      = nghost;
189:   v->petscnative = PETSC_TRUE;

191:   PetscLayoutSetUp(v->map);

193:   s->array           = (PetscScalar*)array;
194:   s->array_allocated = 0;
195:   if (alloc && !array) {
196:     PetscInt n = v->map->n+nghost;
197:     PetscMalloc(n*sizeof(PetscScalar),&s->array);
198:     PetscLogObjectMemory(v,n*sizeof(PetscScalar));
199:     PetscMemzero(s->array,v->map->n*sizeof(PetscScalar));
200:     s->array_allocated = s->array;
201:   }

203:   /* By default parallel vectors do not have local representation */
204:   s->localrep    = 0;
205:   s->localupdate = 0;

207:   v->stash.insertmode = NOT_SET_VALUES;
208:   /* create the stashes. The block-size for bstash is set later when
209:      VecSetValuesBlocked is called.
210:   */
211:   VecStashCreate_Private(PetscObjectComm((PetscObject)v),1,&v->stash);
212:   VecStashCreate_Private(PetscObjectComm((PetscObject)v),v->map->bs,&v->bstash);

214: #if defined(PETSC_HAVE_MATLAB_ENGINE)
215:   PetscObjectComposeFunction((PetscObject)v,"PetscMatlabEnginePut_C",VecMatlabEnginePut_Default);
216:   PetscObjectComposeFunction((PetscObject)v,"PetscMatlabEngineGet_C",VecMatlabEngineGet_Default);
217: #endif
218:   PetscObjectChangeTypeName((PetscObject)v,VECMPI);
219:   return(0);
220: }

222: /*MC
223:    VECMPI - VECMPI = "mpi" - The basic parallel vector

225:    Options Database Keys:
226: . -vec_type mpi - sets the vector type to VECMPI during a call to VecSetFromOptions()

228:   Level: beginner

230: .seealso: VecCreate(), VecSetType(), VecSetFromOptions(), VecCreateMpiWithArray(), VECMPI, VecType, VecCreateMPI(), VecCreateMpi()
231: M*/

235: PETSC_EXTERN PetscErrorCode VecCreate_MPI(Vec vv)
236: {

240:   VecCreate_MPI_Private(vv,PETSC_TRUE,0,0);
241:   return(0);
242: }

244: /*MC
245:    VECSTANDARD = "standard" - A VECSEQ on one process and VECMPI on more than one process

247:    Options Database Keys:
248: . -vec_type standard - sets a vector type to standard on calls to VecSetFromOptions()

250:   Level: beginner

252: .seealso: VecCreateSeq(), VecCreateMPI()
253: M*/

257: PETSC_EXTERN PetscErrorCode VecCreate_Standard(Vec v)
258: {
260:   PetscMPIInt    size;

263:   MPI_Comm_size(PetscObjectComm((PetscObject)v),&size);
264:   if (size == 1) {
265:     VecSetType(v,VECSEQ);
266:   } else {
267:     VecSetType(v,VECMPI);
268:   }
269:   return(0);
270: }

274: /*@C
275:    VecCreateMPIWithArray - Creates a parallel, array-style vector,
276:    where the user provides the array space to store the vector values.

278:    Collective on MPI_Comm

280:    Input Parameters:
281: +  comm  - the MPI communicator to use
282: .  bs    - block size, same meaning as VecSetBlockSize()
283: .  n     - local vector length, cannot be PETSC_DECIDE
284: .  N     - global vector length (or PETSC_DECIDE to have calculated)
285: -  array - the user provided array to store the vector values

287:    Output Parameter:
288: .  vv - the vector

290:    Notes:
291:    Use VecDuplicate() or VecDuplicateVecs() to form additional vectors of the
292:    same type as an existing vector.

294:    If the user-provided array is NULL, then VecPlaceArray() can be used
295:    at a later stage to SET the array for storing the vector values.

297:    PETSc does NOT free the array when the vector is destroyed via VecDestroy().
298:    The user should not free the array until the vector is destroyed.

300:    Level: intermediate

302:    Concepts: vectors^creating with array

304: .seealso: VecCreateSeqWithArray(), VecCreate(), VecDuplicate(), VecDuplicateVecs(), VecCreateGhost(),
305:           VecCreateMPI(), VecCreateGhostWithArray(), VecPlaceArray()

307: @*/
308: PetscErrorCode  VecCreateMPIWithArray(MPI_Comm comm,PetscInt bs,PetscInt n,PetscInt N,const PetscScalar array[],Vec *vv)
309: {

313:   if (n == PETSC_DECIDE) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Must set local size of vector");
314:   PetscSplitOwnership(comm,&n,&N);
315:   VecCreate(comm,vv);
316:   VecSetSizes(*vv,n,N);
317:   VecSetBlockSize(*vv,bs);
318:   VecCreate_MPI_Private(*vv,PETSC_FALSE,0,array);
319:   return(0);
320: }

324: /*@C
325:    VecCreateGhostWithArray - Creates a parallel vector with ghost padding on each processor;
326:    the caller allocates the array space.

328:    Collective on MPI_Comm

330:    Input Parameters:
331: +  comm - the MPI communicator to use
332: .  n - local vector length
333: .  N - global vector length (or PETSC_DECIDE to have calculated if n is given)
334: .  nghost - number of local ghost points
335: .  ghosts - global indices of ghost points (or NULL if not needed), these do not need to be in increasing order (sorted)
336: -  array - the space to store the vector values (as long as n + nghost)

338:    Output Parameter:
339: .  vv - the global vector representation (without ghost points as part of vector)

341:    Notes:
342:    Use VecGhostGetLocalForm() to access the local, ghosted representation
343:    of the vector.

345:    This also automatically sets the ISLocalToGlobalMapping() for this vector.

347:    Level: advanced

349:    Concepts: vectors^creating with array
350:    Concepts: vectors^ghosted

352: .seealso: VecCreate(), VecGhostGetLocalForm(), VecGhostRestoreLocalForm(),
353:           VecCreateGhost(), VecCreateSeqWithArray(), VecCreateMPIWithArray(),
354:           VecCreateGhostBlock(), VecCreateGhostBlockWithArray(), VecMPISetGhost()

356: @*/
357: PetscErrorCode  VecCreateGhostWithArray(MPI_Comm comm,PetscInt n,PetscInt N,PetscInt nghost,const PetscInt ghosts[],const PetscScalar array[],Vec *vv)
358: {
359:   PetscErrorCode         ierr;
360:   Vec_MPI                *w;
361:   PetscScalar            *larray;
362:   IS                     from,to;
363:   ISLocalToGlobalMapping ltog;
364:   PetscInt               rstart,i,*indices;

367:   *vv = 0;

369:   if (n == PETSC_DECIDE)      SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Must set local size");
370:   if (nghost == PETSC_DECIDE) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Must set local ghost size");
371:   if (nghost < 0)             SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Ghost length must be >= 0");
372:   PetscSplitOwnership(comm,&n,&N);
373:   /* Create global representation */
374:   VecCreate(comm,vv);
375:   VecSetSizes(*vv,n,N);
376:   VecCreate_MPI_Private(*vv,PETSC_TRUE,nghost,array);
377:   w    = (Vec_MPI*)(*vv)->data;
378:   /* Create local representation */
379:   VecGetArray(*vv,&larray);
380:   VecCreateSeqWithArray(PETSC_COMM_SELF,1,n+nghost,larray,&w->localrep);
381:   PetscLogObjectParent(*vv,w->localrep);
382:   VecRestoreArray(*vv,&larray);

384:   /*
385:        Create scatter context for scattering (updating) ghost values
386:   */
387:   ISCreateGeneral(comm,nghost,ghosts,PETSC_COPY_VALUES,&from);
388:   ISCreateStride(PETSC_COMM_SELF,nghost,n,1,&to);
389:   VecScatterCreate(*vv,from,w->localrep,to,&w->localupdate);
390:   PetscLogObjectParent(*vv,w->localupdate);
391:   ISDestroy(&to);
392:   ISDestroy(&from);

394:   /* set local to global mapping for ghosted vector */
395:   PetscMalloc((n+nghost)*sizeof(PetscInt),&indices);
396:   VecGetOwnershipRange(*vv,&rstart,NULL);
397:   for (i=0; i<n; i++) {
398:     indices[i] = rstart + i;
399:   }
400:   for (i=0; i<nghost; i++) {
401:     indices[n+i] = ghosts[i];
402:   }
403:   ISLocalToGlobalMappingCreate(comm,n+nghost,indices,PETSC_OWN_POINTER,&ltog);
404:   VecSetLocalToGlobalMapping(*vv,ltog);
405:   ISLocalToGlobalMappingDestroy(&ltog);
406:   return(0);
407: }

411: /*@
412:    VecCreateGhost - Creates a parallel vector with ghost padding on each processor.

414:    Collective on MPI_Comm

416:    Input Parameters:
417: +  comm - the MPI communicator to use
418: .  n - local vector length
419: .  N - global vector length (or PETSC_DECIDE to have calculated if n is given)
420: .  nghost - number of local ghost points
421: -  ghosts - global indices of ghost points, these do not need to be in increasing order (sorted)

423:    Output Parameter:
424: .  vv - the global vector representation (without ghost points as part of vector)

426:    Notes:
427:    Use VecGhostGetLocalForm() to access the local, ghosted representation
428:    of the vector.

430:    This also automatically sets the ISLocalToGlobalMapping() for this vector.

432:    Level: advanced

434:    Concepts: vectors^ghosted

436: .seealso: VecCreateSeq(), VecCreate(), VecDuplicate(), VecDuplicateVecs(), VecCreateMPI(),
437:           VecGhostGetLocalForm(), VecGhostRestoreLocalForm(), VecGhostUpdateBegin(),
438:           VecCreateGhostWithArray(), VecCreateMPIWithArray(), VecGhostUpdateEnd(),
439:           VecCreateGhostBlock(), VecCreateGhostBlockWithArray(), VecMPISetGhost()

441: @*/
442: PetscErrorCode  VecCreateGhost(MPI_Comm comm,PetscInt n,PetscInt N,PetscInt nghost,const PetscInt ghosts[],Vec *vv)
443: {

447:   VecCreateGhostWithArray(comm,n,N,nghost,ghosts,0,vv);
448:   return(0);
449: }

453: /*@
454:    VecMPISetGhost - Sets the ghost points for an MPI ghost vector

456:    Collective on Vec

458:    Input Parameters:
459: +  vv - the MPI vector
460: .  nghost - number of local ghost points
461: -  ghosts - global indices of ghost points, these do not need to be in increasing order (sorted)


464:    Notes:
465:    Use VecGhostGetLocalForm() to access the local, ghosted representation
466:    of the vector.

468:    This also automatically sets the ISLocalToGlobalMapping() for this vector.

470:    You must call this AFTER you have set the type of the vector (with VecSetType()) and the size (with VecSetSizes()).

472:    Level: advanced

474:    Concepts: vectors^ghosted

476: .seealso: VecCreateSeq(), VecCreate(), VecDuplicate(), VecDuplicateVecs(), VecCreateMPI(),
477:           VecGhostGetLocalForm(), VecGhostRestoreLocalForm(), VecGhostUpdateBegin(),
478:           VecCreateGhostWithArray(), VecCreateMPIWithArray(), VecGhostUpdateEnd(),
479:           VecCreateGhostBlock(), VecCreateGhostBlockWithArray()

481: @*/
482: PetscErrorCode  VecMPISetGhost(Vec vv,PetscInt nghost,const PetscInt ghosts[])
483: {
485:   PetscBool      flg;

488:   PetscObjectTypeCompare((PetscObject)vv,VECMPI,&flg);
489:   /* if already fully existant VECMPI then basically destroy it and rebuild with ghosting */
490:   if (flg) {
491:     PetscInt               n,N;
492:     Vec_MPI                *w;
493:     PetscScalar            *larray;
494:     IS                     from,to;
495:     ISLocalToGlobalMapping ltog;
496:     PetscInt               rstart,i,*indices;
497:     MPI_Comm               comm;

499:     PetscObjectGetComm((PetscObject)vv,&comm);
500:     n    = vv->map->n;
501:     N    = vv->map->N;
502:     (*vv->ops->destroy)(vv);
503:     VecSetSizes(vv,n,N);
504:     VecCreate_MPI_Private(vv,PETSC_TRUE,nghost,NULL);
505:     w    = (Vec_MPI*)(vv)->data;
506:     /* Create local representation */
507:     VecGetArray(vv,&larray);
508:     VecCreateSeqWithArray(PETSC_COMM_SELF,1,n+nghost,larray,&w->localrep);
509:     PetscLogObjectParent(vv,w->localrep);
510:     VecRestoreArray(vv,&larray);

512:     /*
513:      Create scatter context for scattering (updating) ghost values
514:      */
515:     ISCreateGeneral(comm,nghost,ghosts,PETSC_COPY_VALUES,&from);
516:     ISCreateStride(PETSC_COMM_SELF,nghost,n,1,&to);
517:     VecScatterCreate(vv,from,w->localrep,to,&w->localupdate);
518:     PetscLogObjectParent(vv,w->localupdate);
519:     ISDestroy(&to);
520:     ISDestroy(&from);

522:     /* set local to global mapping for ghosted vector */
523:     PetscMalloc((n+nghost)*sizeof(PetscInt),&indices);
524:     VecGetOwnershipRange(vv,&rstart,NULL);

526:     for (i=0; i<n; i++)      indices[i]   = rstart + i;
527:     for (i=0; i<nghost; i++) indices[n+i] = ghosts[i];

529:     ISLocalToGlobalMappingCreate(comm,n+nghost,indices,PETSC_OWN_POINTER,&ltog);
530:     VecSetLocalToGlobalMapping(vv,ltog);
531:     ISLocalToGlobalMappingDestroy(&ltog);
532:   } else if (vv->ops->create == VecCreate_MPI) SETERRQ(PetscObjectComm((PetscObject)vv),PETSC_ERR_ARG_WRONGSTATE,"Must set local or global size before setting ghosting");
533:   else if (!((PetscObject)vv)->type_name) SETERRQ(PetscObjectComm((PetscObject)vv),PETSC_ERR_ARG_WRONGSTATE,"Must set type to VECMPI before ghosting");
534:   return(0);
535: }


538: /* ------------------------------------------------------------------------------------------*/
541: /*@C
542:    VecCreateGhostBlockWithArray - Creates a parallel vector with ghost padding on each processor;
543:    the caller allocates the array space. Indices in the ghost region are based on blocks.

545:    Collective on MPI_Comm

547:    Input Parameters:
548: +  comm - the MPI communicator to use
549: .  bs - block size
550: .  n - local vector length
551: .  N - global vector length (or PETSC_DECIDE to have calculated if n is given)
552: .  nghost - number of local ghost blocks
553: .  ghosts - global indices of ghost blocks (or NULL if not needed), counts are by block not by index, these do not need to be in increasing order (sorted)
554: -  array - the space to store the vector values (as long as n + nghost*bs)

556:    Output Parameter:
557: .  vv - the global vector representation (without ghost points as part of vector)

559:    Notes:
560:    Use VecGhostGetLocalForm() to access the local, ghosted representation
561:    of the vector.

563:    n is the local vector size (total local size not the number of blocks) while nghost
564:    is the number of blocks in the ghost portion, i.e. the number of elements in the ghost
565:    portion is bs*nghost

567:    Level: advanced

569:    Concepts: vectors^creating ghosted
570:    Concepts: vectors^creating with array

572: .seealso: VecCreate(), VecGhostGetLocalForm(), VecGhostRestoreLocalForm(),
573:           VecCreateGhost(), VecCreateSeqWithArray(), VecCreateMPIWithArray(),
574:           VecCreateGhostWithArray(), VecCreateGhostBlock()

576: @*/
577: PetscErrorCode  VecCreateGhostBlockWithArray(MPI_Comm comm,PetscInt bs,PetscInt n,PetscInt N,PetscInt nghost,const PetscInt ghosts[],const PetscScalar array[],Vec *vv)
578: {
579:   PetscErrorCode         ierr;
580:   Vec_MPI                *w;
581:   PetscScalar            *larray;
582:   IS                     from,to;
583:   ISLocalToGlobalMapping ltog;
584:   PetscInt               rstart,i,nb,*indices;

587:   *vv = 0;

589:   if (n == PETSC_DECIDE)      SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Must set local size");
590:   if (nghost == PETSC_DECIDE) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Must set local ghost size");
591:   if (nghost < 0)             SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Ghost length must be >= 0");
592:   if (n % bs)                 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Local size must be a multiple of block size");
593:   PetscSplitOwnership(comm,&n,&N);
594:   /* Create global representation */
595:   VecCreate(comm,vv);
596:   VecSetSizes(*vv,n,N);
597:   VecSetBlockSize(*vv,bs);
598:   VecCreate_MPI_Private(*vv,PETSC_TRUE,nghost*bs,array);
599:   w    = (Vec_MPI*)(*vv)->data;
600:   /* Create local representation */
601:   VecGetArray(*vv,&larray);
602:   VecCreateSeqWithArray(PETSC_COMM_SELF,bs,n+bs*nghost,larray,&w->localrep);
603:   PetscLogObjectParent(*vv,w->localrep);
604:   VecRestoreArray(*vv,&larray);

606:   /*
607:        Create scatter context for scattering (updating) ghost values
608:   */
609:   ISCreateBlock(comm,bs,nghost,ghosts,PETSC_COPY_VALUES,&from);
610:   ISCreateStride(PETSC_COMM_SELF,bs*nghost,n,1,&to);
611:   VecScatterCreate(*vv,from,w->localrep,to,&w->localupdate);
612:   PetscLogObjectParent(*vv,w->localupdate);
613:   ISDestroy(&to);
614:   ISDestroy(&from);

616:   /* set local to global mapping for ghosted vector */
617:   nb   = n/bs;
618:   PetscMalloc((nb+nghost)*sizeof(PetscInt),&indices);
619:   VecGetOwnershipRange(*vv,&rstart,NULL);

621:   for (i=0; i<nb; i++)      indices[i]    = rstart + i*bs;
622:   for (i=0; i<nghost; i++)  indices[nb+i] = ghosts[i];

624:   ISLocalToGlobalMappingCreate(comm,nb+nghost,indices,PETSC_OWN_POINTER,&ltog);
625:   VecSetLocalToGlobalMappingBlock(*vv,ltog);
626:   ISLocalToGlobalMappingDestroy(&ltog);
627:   return(0);
628: }

632: /*@
633:    VecCreateGhostBlock - Creates a parallel vector with ghost padding on each processor.
634:         The indicing of the ghost points is done with blocks.

636:    Collective on MPI_Comm

638:    Input Parameters:
639: +  comm - the MPI communicator to use
640: .  bs - the block size
641: .  n - local vector length
642: .  N - global vector length (or PETSC_DECIDE to have calculated if n is given)
643: .  nghost - number of local ghost blocks
644: -  ghosts - global indices of ghost blocks, counts are by block, not by individual index, these do not need to be in increasing order (sorted)

646:    Output Parameter:
647: .  vv - the global vector representation (without ghost points as part of vector)

649:    Notes:
650:    Use VecGhostGetLocalForm() to access the local, ghosted representation
651:    of the vector.

653:    n is the local vector size (total local size not the number of blocks) while nghost
654:    is the number of blocks in the ghost portion, i.e. the number of elements in the ghost
655:    portion is bs*nghost

657:    Level: advanced

659:    Concepts: vectors^ghosted

661: .seealso: VecCreateSeq(), VecCreate(), VecDuplicate(), VecDuplicateVecs(), VecCreateMPI(),
662:           VecGhostGetLocalForm(), VecGhostRestoreLocalForm(),
663:           VecCreateGhostWithArray(), VecCreateMPIWithArray(), VecCreateGhostBlockWithArray()

665: @*/
666: PetscErrorCode  VecCreateGhostBlock(MPI_Comm comm,PetscInt bs,PetscInt n,PetscInt N,PetscInt nghost,const PetscInt ghosts[],Vec *vv)
667: {

671:   VecCreateGhostBlockWithArray(comm,bs,n,N,nghost,ghosts,0,vv);
672:   return(0);
673: }