Actual source code: pbvec.c

petsc-3.5.2 2014-09-08
Report Typos and Errors
  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,PetscAbs(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((PetscObject)*v,(PetscObject)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   = PetscAbs(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:                                 VecStrideSubSetGather_Default,
168:                                 VecStrideSubSetScatter_Default,
169:                                 0,
170:                                 0
171: };

175: /*
176:     VecCreate_MPI_Private - Basic create routine called by VecCreate_MPI() (i.e. VecCreateMPI()),
177:     VecCreateMPIWithArray(), VecCreate_Shared() (i.e. VecCreateShared()), VecCreateGhost(),
178:     VecDuplicate_MPI(), VecCreateGhostWithArray(), VecDuplicate_MPI(), and VecDuplicate_Shared()

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

189:   PetscNewLog(v,&s);
190:   v->data        = (void*)s;
191:   PetscMemcpy(v->ops,&DvOps,sizeof(DvOps));
192:   s->nghost      = nghost;
193:   v->petscnative = PETSC_TRUE;

195:   PetscLayoutSetUp(v->map);

197:   s->array           = (PetscScalar*)array;
198:   s->array_allocated = 0;
199:   if (alloc && !array) {
200:     PetscInt n = v->map->n+nghost;
201:     PetscMalloc1(n,&s->array);
202:     PetscLogObjectMemory((PetscObject)v,n*sizeof(PetscScalar));
203:     PetscMemzero(s->array,v->map->n*sizeof(PetscScalar));
204:     s->array_allocated = s->array;
205:   }

207:   /* By default parallel vectors do not have local representation */
208:   s->localrep    = 0;
209:   s->localupdate = 0;

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

218: #if defined(PETSC_HAVE_MATLAB_ENGINE)
219:   PetscObjectComposeFunction((PetscObject)v,"PetscMatlabEnginePut_C",VecMatlabEnginePut_Default);
220:   PetscObjectComposeFunction((PetscObject)v,"PetscMatlabEngineGet_C",VecMatlabEngineGet_Default);
221: #endif
222:   PetscObjectChangeTypeName((PetscObject)v,VECMPI);
223:   return(0);
224: }

226: /*MC
227:    VECMPI - VECMPI = "mpi" - The basic parallel vector

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

232:   Level: beginner

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

239: PETSC_EXTERN PetscErrorCode VecCreate_MPI(Vec vv)
240: {

244:   VecCreate_MPI_Private(vv,PETSC_TRUE,0,0);
245:   return(0);
246: }

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

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

254:   Level: beginner

256: .seealso: VecCreateSeq(), VecCreateMPI()
257: M*/

261: PETSC_EXTERN PetscErrorCode VecCreate_Standard(Vec v)
262: {
264:   PetscMPIInt    size;

267:   MPI_Comm_size(PetscObjectComm((PetscObject)v),&size);
268:   if (size == 1) {
269:     VecSetType(v,VECSEQ);
270:   } else {
271:     VecSetType(v,VECMPI);
272:   }
273:   return(0);
274: }

278: /*@C
279:    VecCreateMPIWithArray - Creates a parallel, array-style vector,
280:    where the user provides the array space to store the vector values.

282:    Collective on MPI_Comm

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

291:    Output Parameter:
292: .  vv - the vector

294:    Notes:
295:    Use VecDuplicate() or VecDuplicateVecs() to form additional vectors of the
296:    same type as an existing vector.

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

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

304:    Level: intermediate

306:    Concepts: vectors^creating with array

308: .seealso: VecCreateSeqWithArray(), VecCreate(), VecDuplicate(), VecDuplicateVecs(), VecCreateGhost(),
309:           VecCreateMPI(), VecCreateGhostWithArray(), VecPlaceArray()

311: @*/
312: PetscErrorCode  VecCreateMPIWithArray(MPI_Comm comm,PetscInt bs,PetscInt n,PetscInt N,const PetscScalar array[],Vec *vv)
313: {

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

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

332:    Collective on MPI_Comm

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

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

345:    Notes:
346:    Use VecGhostGetLocalForm() to access the local, ghosted representation
347:    of the vector.

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

351:    Level: advanced

353:    Concepts: vectors^creating with array
354:    Concepts: vectors^ghosted

356: .seealso: VecCreate(), VecGhostGetLocalForm(), VecGhostRestoreLocalForm(),
357:           VecCreateGhost(), VecCreateSeqWithArray(), VecCreateMPIWithArray(),
358:           VecCreateGhostBlock(), VecCreateGhostBlockWithArray(), VecMPISetGhost()

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

371:   *vv = 0;

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

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

398:   /* set local to global mapping for ghosted vector */
399:   PetscMalloc1((n+nghost),&indices);
400:   VecGetOwnershipRange(*vv,&rstart,NULL);
401:   for (i=0; i<n; i++) {
402:     indices[i] = rstart + i;
403:   }
404:   for (i=0; i<nghost; i++) {
405:     indices[n+i] = ghosts[i];
406:   }
407:   ISLocalToGlobalMappingCreate(comm,1,n+nghost,indices,PETSC_OWN_POINTER,&ltog);
408:   VecSetLocalToGlobalMapping(*vv,ltog);
409:   ISLocalToGlobalMappingDestroy(&ltog);
410:   return(0);
411: }

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

418:    Collective on MPI_Comm

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

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

430:    Notes:
431:    Use VecGhostGetLocalForm() to access the local, ghosted representation
432:    of the vector.

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

436:    Level: advanced

438:    Concepts: vectors^ghosted

440: .seealso: VecCreateSeq(), VecCreate(), VecDuplicate(), VecDuplicateVecs(), VecCreateMPI(),
441:           VecGhostGetLocalForm(), VecGhostRestoreLocalForm(), VecGhostUpdateBegin(),
442:           VecCreateGhostWithArray(), VecCreateMPIWithArray(), VecGhostUpdateEnd(),
443:           VecCreateGhostBlock(), VecCreateGhostBlockWithArray(), VecMPISetGhost()

445: @*/
446: PetscErrorCode  VecCreateGhost(MPI_Comm comm,PetscInt n,PetscInt N,PetscInt nghost,const PetscInt ghosts[],Vec *vv)
447: {

451:   VecCreateGhostWithArray(comm,n,N,nghost,ghosts,0,vv);
452:   return(0);
453: }

457: /*@
458:    VecMPISetGhost - Sets the ghost points for an MPI ghost vector

460:    Collective on Vec

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


468:    Notes:
469:    Use VecGhostGetLocalForm() to access the local, ghosted representation
470:    of the vector.

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

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

476:    Level: advanced

478:    Concepts: vectors^ghosted

480: .seealso: VecCreateSeq(), VecCreate(), VecDuplicate(), VecDuplicateVecs(), VecCreateMPI(),
481:           VecGhostGetLocalForm(), VecGhostRestoreLocalForm(), VecGhostUpdateBegin(),
482:           VecCreateGhostWithArray(), VecCreateMPIWithArray(), VecGhostUpdateEnd(),
483:           VecCreateGhostBlock(), VecCreateGhostBlockWithArray()

485: @*/
486: PetscErrorCode  VecMPISetGhost(Vec vv,PetscInt nghost,const PetscInt ghosts[])
487: {
489:   PetscBool      flg;

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

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

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

526:     /* set local to global mapping for ghosted vector */
527:     PetscMalloc1((n+nghost),&indices);
528:     VecGetOwnershipRange(vv,&rstart,NULL);

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

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


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

549:    Collective on MPI_Comm

551:    Input Parameters:
552: +  comm - the MPI communicator to use
553: .  bs - block size
554: .  n - local vector length
555: .  N - global vector length (or PETSC_DECIDE to have calculated if n is given)
556: .  nghost - number of local ghost blocks
557: .  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)
558: -  array - the space to store the vector values (as long as n + nghost*bs)

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

563:    Notes:
564:    Use VecGhostGetLocalForm() to access the local, ghosted representation
565:    of the vector.

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

571:    Level: advanced

573:    Concepts: vectors^creating ghosted
574:    Concepts: vectors^creating with array

576: .seealso: VecCreate(), VecGhostGetLocalForm(), VecGhostRestoreLocalForm(),
577:           VecCreateGhost(), VecCreateSeqWithArray(), VecCreateMPIWithArray(),
578:           VecCreateGhostWithArray(), VecCreateGhostBlock()

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

591:   *vv = 0;

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

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

620:   /* set local to global mapping for ghosted vector */
621:   nb   = n/bs;
622:   PetscMalloc1((nb+nghost),&indices);
623:   VecGetOwnershipRange(*vv,&rstart,NULL);

625:   for (i=0; i<nb; i++)      indices[i]    = rstart + i;
626:   for (i=0; i<nghost; i++)  indices[nb+i] = ghosts[i];

628:   ISLocalToGlobalMappingCreate(comm,bs,nb+nghost,indices,PETSC_OWN_POINTER,&ltog);
629:   VecSetLocalToGlobalMapping(*vv,ltog);
630:   ISLocalToGlobalMappingDestroy(&ltog);
631:   return(0);
632: }

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

640:    Collective on MPI_Comm

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

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

653:    Notes:
654:    Use VecGhostGetLocalForm() to access the local, ghosted representation
655:    of the vector.

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

661:    Level: advanced

663:    Concepts: vectors^ghosted

665: .seealso: VecCreateSeq(), VecCreate(), VecDuplicate(), VecDuplicateVecs(), VecCreateMPI(),
666:           VecGhostGetLocalForm(), VecGhostRestoreLocalForm(),
667:           VecCreateGhostWithArray(), VecCreateMPIWithArray(), VecCreateGhostBlockWithArray()

669: @*/
670: PetscErrorCode  VecCreateGhostBlock(MPI_Comm comm,PetscInt bs,PetscInt n,PetscInt N,PetscInt nghost,const PetscInt ghosts[],Vec *vv)
671: {

675:   VecCreateGhostBlockWithArray(comm,bs,n,N,nghost,ghosts,0,vv);
676:   return(0);
677: }