1 2 /* 3 This file contains routines for Parallel vector operations. 4 */ 5 #include <petscoptions.h> 6 #include <../src/vec/vec/impls/mpi/pvecimpl.h> /*I "petscvec.h" I*/ 7 8 #undef __FUNCT__ 9 #define __FUNCT__ "VecDot_MPI" 10 static PetscErrorCode VecDot_MPI(Vec xin,Vec yin,PetscScalar *z) 11 { 12 PetscScalar sum,work; 13 PetscErrorCode ierr; 14 15 PetscFunctionBegin; 16 ierr = VecDot_Seq(xin,yin,&work);CHKERRQ(ierr); 17 ierr = MPIU_Allreduce(&work,&sum,1,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)xin));CHKERRQ(ierr); 18 *z = sum; 19 PetscFunctionReturn(0); 20 } 21 22 #undef __FUNCT__ 23 #define __FUNCT__ "VecTDot_MPI" 24 static PetscErrorCode VecTDot_MPI(Vec xin,Vec yin,PetscScalar *z) 25 { 26 PetscScalar sum,work; 27 PetscErrorCode ierr; 28 29 PetscFunctionBegin; 30 ierr = VecTDot_Seq(xin,yin,&work);CHKERRQ(ierr); 31 ierr = MPIU_Allreduce(&work,&sum,1,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)xin));CHKERRQ(ierr); 32 *z = sum; 33 PetscFunctionReturn(0); 34 } 35 36 extern PetscErrorCode VecView_MPI_Draw(Vec,PetscViewer); 37 38 #undef __FUNCT__ 39 #define __FUNCT__ "VecPlaceArray_MPI" 40 static PetscErrorCode VecPlaceArray_MPI(Vec vin,const PetscScalar *a) 41 { 42 PetscErrorCode ierr; 43 Vec_MPI *v = (Vec_MPI*)vin->data; 44 45 PetscFunctionBegin; 46 if (v->unplacedarray) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"VecPlaceArray() was already called on this vector, without a call to VecResetArray()"); 47 v->unplacedarray = v->array; /* save previous array so reset can bring it back */ 48 v->array = (PetscScalar*)a; 49 if (v->localrep) { 50 ierr = VecPlaceArray(v->localrep,a);CHKERRQ(ierr); 51 } 52 PetscFunctionReturn(0); 53 } 54 55 extern PetscErrorCode VecGetValues_MPI(Vec,PetscInt,const PetscInt[],PetscScalar[]); 56 57 #undef __FUNCT__ 58 #define __FUNCT__ "VecDuplicate_MPI" 59 static PetscErrorCode VecDuplicate_MPI(Vec win,Vec *v) 60 { 61 PetscErrorCode ierr; 62 Vec_MPI *vw,*w = (Vec_MPI*)win->data; 63 PetscScalar *array; 64 65 PetscFunctionBegin; 66 ierr = VecCreate(PetscObjectComm((PetscObject)win),v);CHKERRQ(ierr); 67 ierr = PetscLayoutReference(win->map,&(*v)->map);CHKERRQ(ierr); 68 69 ierr = VecCreate_MPI_Private(*v,PETSC_TRUE,w->nghost,0);CHKERRQ(ierr); 70 vw = (Vec_MPI*)(*v)->data; 71 ierr = PetscMemcpy((*v)->ops,win->ops,sizeof(struct _VecOps));CHKERRQ(ierr); 72 73 /* save local representation of the parallel vector (and scatter) if it exists */ 74 if (w->localrep) { 75 ierr = VecGetArray(*v,&array);CHKERRQ(ierr); 76 ierr = VecCreateSeqWithArray(PETSC_COMM_SELF,PetscAbs(win->map->bs),win->map->n+w->nghost,array,&vw->localrep);CHKERRQ(ierr); 77 ierr = PetscMemcpy(vw->localrep->ops,w->localrep->ops,sizeof(struct _VecOps));CHKERRQ(ierr); 78 ierr = VecRestoreArray(*v,&array);CHKERRQ(ierr); 79 ierr = PetscLogObjectParent((PetscObject)*v,(PetscObject)vw->localrep);CHKERRQ(ierr); 80 81 vw->localupdate = w->localupdate; 82 if (vw->localupdate) { 83 ierr = PetscObjectReference((PetscObject)vw->localupdate);CHKERRQ(ierr); 84 } 85 } 86 87 /* New vector should inherit stashing property of parent */ 88 (*v)->stash.donotstash = win->stash.donotstash; 89 (*v)->stash.ignorenegidx = win->stash.ignorenegidx; 90 91 ierr = PetscObjectListDuplicate(((PetscObject)win)->olist,&((PetscObject)(*v))->olist);CHKERRQ(ierr); 92 ierr = PetscFunctionListDuplicate(((PetscObject)win)->qlist,&((PetscObject)(*v))->qlist);CHKERRQ(ierr); 93 94 (*v)->map->bs = PetscAbs(win->map->bs); 95 (*v)->bstash.bs = win->bstash.bs; 96 PetscFunctionReturn(0); 97 } 98 99 100 #undef __FUNCT__ 101 #define __FUNCT__ "VecSetOption_MPI" 102 static PetscErrorCode VecSetOption_MPI(Vec V,VecOption op,PetscBool flag) 103 { 104 Vec_MPI *v = (Vec_MPI*)V->data; 105 PetscFunctionBegin; 106 switch (op) { 107 case VEC_IGNORE_OFF_PROC_ENTRIES: V->stash.donotstash = flag; 108 break; 109 case VEC_IGNORE_NEGATIVE_INDICES: V->stash.ignorenegidx = flag; 110 break; 111 case VEC_SUBSET_OFF_PROC_ENTRIES: v->assembly_subset = flag; 112 break; 113 } 114 PetscFunctionReturn(0); 115 } 116 117 118 #undef __FUNCT__ 119 #define __FUNCT__ "VecResetArray_MPI" 120 static PetscErrorCode VecResetArray_MPI(Vec vin) 121 { 122 Vec_MPI *v = (Vec_MPI*)vin->data; 123 PetscErrorCode ierr; 124 125 PetscFunctionBegin; 126 v->array = v->unplacedarray; 127 v->unplacedarray = 0; 128 if (v->localrep) { 129 ierr = VecResetArray(v->localrep);CHKERRQ(ierr); 130 } 131 PetscFunctionReturn(0); 132 } 133 134 #undef __FUNCT__ 135 #define __FUNCT__ "VecAssemblySend_MPI_Private" 136 static PetscErrorCode VecAssemblySend_MPI_Private(MPI_Comm comm,const PetscMPIInt tag[],PetscMPIInt rankid,PetscMPIInt rank,void *sdata,MPI_Request req[],void *ctx) 137 { 138 Vec X = (Vec)ctx; 139 Vec_MPI *x = (Vec_MPI*)X->data; 140 VecAssemblyHeader *hdr = (VecAssemblyHeader*)sdata; 141 PetscInt bs = X->map->bs; 142 PetscErrorCode ierr; 143 144 PetscFunctionBegin; 145 /* x->recvhdr only exists when we are reusing a communication network. In that case, some messages can be empty, but 146 * we have to send them this time if we sent them before because the receiver is expecting them. */ 147 if (hdr->count || (x->recvhdr && x->sendptrs[rankid].ints)) { 148 ierr = MPI_Isend(x->sendptrs[rankid].ints,hdr->count,MPIU_INT,rank,tag[0],comm,&req[0]);CHKERRQ(ierr); 149 ierr = MPI_Isend(x->sendptrs[rankid].scalars,hdr->count,MPIU_SCALAR,rank,tag[1],comm,&req[1]);CHKERRQ(ierr); 150 } 151 if (hdr->bcount || (x->recvhdr && x->sendptrs[rankid].intb)) { 152 ierr = MPI_Isend(x->sendptrs[rankid].intb,hdr->bcount,MPIU_INT,rank,tag[2],comm,&req[2]);CHKERRQ(ierr); 153 ierr = MPI_Isend(x->sendptrs[rankid].scalarb,hdr->bcount*bs,MPIU_SCALAR,rank,tag[3],comm,&req[3]);CHKERRQ(ierr); 154 } 155 PetscFunctionReturn(0); 156 } 157 158 #undef __FUNCT__ 159 #define __FUNCT__ "VecAssemblyRecv_MPI_Private" 160 static PetscErrorCode VecAssemblyRecv_MPI_Private(MPI_Comm comm,const PetscMPIInt tag[],PetscMPIInt rank,void *rdata,MPI_Request req[],void *ctx) 161 { 162 Vec X = (Vec)ctx; 163 Vec_MPI *x = (Vec_MPI*)X->data; 164 VecAssemblyHeader *hdr = (VecAssemblyHeader*)rdata; 165 PetscErrorCode ierr; 166 PetscInt bs = X->map->bs; 167 VecAssemblyFrame *frame; 168 169 PetscFunctionBegin; 170 ierr = PetscSegBufferGet(x->segrecvframe,1,&frame);CHKERRQ(ierr); 171 172 if (hdr->count) { 173 ierr = PetscSegBufferGet(x->segrecvint,hdr->count,&frame->ints);CHKERRQ(ierr); 174 ierr = MPI_Irecv(frame->ints,hdr->count,MPIU_INT,rank,tag[0],comm,&req[0]);CHKERRQ(ierr); 175 ierr = PetscSegBufferGet(x->segrecvscalar,hdr->count,&frame->scalars);CHKERRQ(ierr); 176 ierr = MPI_Irecv(frame->scalars,hdr->count,MPIU_SCALAR,rank,tag[1],comm,&req[1]);CHKERRQ(ierr); 177 frame->pendings = 2; 178 } else { 179 frame->ints = NULL; 180 frame->scalars = NULL; 181 frame->pendings = 0; 182 } 183 184 if (hdr->bcount) { 185 ierr = PetscSegBufferGet(x->segrecvint,hdr->bcount,&frame->intb);CHKERRQ(ierr); 186 ierr = MPI_Irecv(frame->intb,hdr->bcount,MPIU_INT,rank,tag[2],comm,&req[2]);CHKERRQ(ierr); 187 ierr = PetscSegBufferGet(x->segrecvscalar,hdr->bcount*bs,&frame->scalarb);CHKERRQ(ierr); 188 ierr = MPI_Irecv(frame->scalarb,hdr->bcount*bs,MPIU_SCALAR,rank,tag[3],comm,&req[3]);CHKERRQ(ierr); 189 frame->pendingb = 2; 190 } else { 191 frame->intb = NULL; 192 frame->scalarb = NULL; 193 frame->pendingb = 0; 194 } 195 PetscFunctionReturn(0); 196 } 197 198 #undef __FUNCT__ 199 #define __FUNCT__ "VecAssemblyBegin_MPI_BTS" 200 static PetscErrorCode VecAssemblyBegin_MPI_BTS(Vec X) 201 { 202 Vec_MPI *x = (Vec_MPI*)X->data; 203 PetscErrorCode ierr; 204 MPI_Comm comm; 205 PetscInt i,j,jb,bs; 206 207 PetscFunctionBegin; 208 if (X->stash.donotstash) PetscFunctionReturn(0); 209 210 ierr = PetscObjectGetComm((PetscObject)X,&comm);CHKERRQ(ierr); 211 ierr = VecGetBlockSize(X,&bs);CHKERRQ(ierr); 212 #if defined(PETSC_USE_DEBUG) 213 { 214 InsertMode addv; 215 ierr = MPIU_Allreduce((PetscEnum*)&X->stash.insertmode,(PetscEnum*)&addv,1,MPIU_ENUM,MPI_BOR,comm);CHKERRQ(ierr); 216 if (addv == (ADD_VALUES|INSERT_VALUES)) SETERRQ(comm,PETSC_ERR_ARG_NOTSAMETYPE,"Some processors inserted values while others added"); 217 } 218 #endif 219 X->bstash.insertmode = X->stash.insertmode; /* Block stash implicitly tracks InsertMode of scalar stash */ 220 221 ierr = VecStashSortCompress_Private(&X->stash);CHKERRQ(ierr); 222 ierr = VecStashSortCompress_Private(&X->bstash);CHKERRQ(ierr); 223 224 if (!x->sendranks) { 225 PetscMPIInt nowners,bnowners,*owners,*bowners; 226 PetscInt ntmp; 227 ierr = VecStashGetOwnerList_Private(&X->stash,X->map,&nowners,&owners);CHKERRQ(ierr); 228 ierr = VecStashGetOwnerList_Private(&X->bstash,X->map,&bnowners,&bowners);CHKERRQ(ierr); 229 ierr = PetscMergeMPIIntArray(nowners,owners,bnowners,bowners,&ntmp,&x->sendranks);CHKERRQ(ierr); 230 x->nsendranks = ntmp; 231 ierr = PetscFree(owners);CHKERRQ(ierr); 232 ierr = PetscFree(bowners);CHKERRQ(ierr); 233 ierr = PetscMalloc1(x->nsendranks,&x->sendhdr);CHKERRQ(ierr); 234 ierr = PetscCalloc1(x->nsendranks,&x->sendptrs);CHKERRQ(ierr); 235 } 236 for (i=0,j=0,jb=0; i<x->nsendranks; i++) { 237 PetscMPIInt rank = x->sendranks[i]; 238 x->sendhdr[i].insertmode = X->stash.insertmode; 239 /* Initialize pointers for non-empty stashes the first time around. Subsequent assemblies with 240 * VEC_SUBSET_OFF_PROC_ENTRIES will leave the old pointers (dangling because the stash has been collected) when 241 * there is nothing new to send, so that size-zero messages get sent instead. */ 242 x->sendhdr[i].count = 0; 243 if (X->stash.n) { 244 x->sendptrs[i].ints = &X->stash.idx[j]; 245 x->sendptrs[i].scalars = &X->stash.array[j]; 246 for ( ; j<X->stash.n && X->stash.idx[j] < X->map->range[rank+1]; j++) x->sendhdr[i].count++; 247 } 248 x->sendhdr[i].bcount = 0; 249 if (X->bstash.n) { 250 x->sendptrs[i].intb = &X->bstash.idx[jb]; 251 x->sendptrs[i].scalarb = &X->bstash.array[jb*bs]; 252 for ( ; jb<X->bstash.n && X->bstash.idx[jb]*bs < X->map->range[rank+1]; jb++) x->sendhdr[i].bcount++; 253 } 254 } 255 256 if (!x->segrecvint) {ierr = PetscSegBufferCreate(sizeof(PetscInt),1000,&x->segrecvint);CHKERRQ(ierr);} 257 if (!x->segrecvscalar) {ierr = PetscSegBufferCreate(sizeof(PetscScalar),1000,&x->segrecvscalar);CHKERRQ(ierr);} 258 if (!x->segrecvframe) {ierr = PetscSegBufferCreate(sizeof(VecAssemblyFrame),50,&x->segrecvframe);CHKERRQ(ierr);} 259 if (x->recvhdr) { /* VEC_SUBSET_OFF_PROC_ENTRIES and this is not the first assembly */ 260 PetscMPIInt tag[4]; 261 if (!x->assembly_subset) SETERRQ(comm,PETSC_ERR_PLIB,"Attempt to reuse rendezvous when not VEC_SUBSET_OFF_PROC_ENTRIES"); 262 for (i=0; i<4; i++) {ierr = PetscCommGetNewTag(comm,&tag[i]);CHKERRQ(ierr);} 263 for (i=0; i<x->nsendranks; i++) { 264 ierr = VecAssemblySend_MPI_Private(comm,tag,i,x->sendranks[i],x->sendhdr+i,x->sendreqs+4*i,X);CHKERRQ(ierr); 265 } 266 for (i=0; i<x->nrecvranks; i++) { 267 ierr = VecAssemblyRecv_MPI_Private(comm,tag,x->recvranks[i],x->recvhdr+i,x->recvreqs+4*i,X);CHKERRQ(ierr); 268 } 269 x->use_status = PETSC_TRUE; 270 } else { /* First time */ 271 ierr = PetscCommBuildTwoSidedFReq(comm,3,MPIU_INT,x->nsendranks,x->sendranks,(PetscInt*)x->sendhdr,&x->nrecvranks,&x->recvranks,&x->recvhdr,4,&x->sendreqs,&x->recvreqs,VecAssemblySend_MPI_Private,VecAssemblyRecv_MPI_Private,X);CHKERRQ(ierr); 272 x->use_status = PETSC_FALSE; 273 } 274 275 { 276 PetscInt nstash,reallocs; 277 ierr = VecStashGetInfo_Private(&X->stash,&nstash,&reallocs);CHKERRQ(ierr); 278 ierr = PetscInfo2(X,"Stash has %D entries, uses %D mallocs.\n",nstash,reallocs);CHKERRQ(ierr); 279 ierr = VecStashGetInfo_Private(&X->bstash,&nstash,&reallocs);CHKERRQ(ierr); 280 ierr = PetscInfo2(X,"Block-Stash has %D entries, uses %D mallocs.\n",nstash,reallocs);CHKERRQ(ierr); 281 } 282 PetscFunctionReturn(0); 283 } 284 285 #undef __FUNCT__ 286 #define __FUNCT__ "VecAssemblyEnd_MPI_BTS" 287 static PetscErrorCode VecAssemblyEnd_MPI_BTS(Vec X) 288 { 289 Vec_MPI *x = (Vec_MPI*)X->data; 290 PetscInt bs = X->map->bs; 291 PetscMPIInt npending,*some_indices,r; 292 MPI_Status *some_statuses; 293 PetscScalar *xarray; 294 PetscErrorCode ierr; 295 VecAssemblyFrame *frame; 296 297 PetscFunctionBegin; 298 if (X->stash.donotstash) { 299 X->stash.insertmode = NOT_SET_VALUES; 300 X->bstash.insertmode = NOT_SET_VALUES; 301 PetscFunctionReturn(0); 302 } 303 304 ierr = VecGetArray(X,&xarray);CHKERRQ(ierr); 305 ierr = PetscSegBufferExtractInPlace(x->segrecvframe,&frame);CHKERRQ(ierr); 306 ierr = PetscMalloc2(4*x->nrecvranks,&some_indices,x->use_status?4*x->nrecvranks:0,&some_statuses);CHKERRQ(ierr); 307 for (r=0,npending=0; r<x->nrecvranks; r++) npending += frame[r].pendings + frame[r].pendingb; 308 while (npending>0) { 309 PetscMPIInt ndone,ii; 310 /* Filling MPI_Status fields requires some resources from the MPI library. We skip it on the first assembly, or 311 * when VEC_SUBSET_OFF_PROC_ENTRIES has not been set, because we could exchange exact sizes in the initial 312 * rendezvous. When the rendezvous is elided, however, we use MPI_Status to get actual message lengths, so that 313 * subsequent assembly can set a proper subset of the values. */ 314 ierr = MPI_Waitsome(4*x->nrecvranks,x->recvreqs,&ndone,some_indices,x->use_status?some_statuses:MPI_STATUSES_IGNORE);CHKERRQ(ierr); 315 for (ii=0; ii<ndone; ii++) { 316 PetscInt i = some_indices[ii]/4,j,k; 317 InsertMode imode = (InsertMode)x->recvhdr[i].insertmode; 318 PetscInt *recvint; 319 PetscScalar *recvscalar; 320 PetscBool intmsg = (PetscBool)(some_indices[ii]%2 == 0); 321 PetscBool blockmsg = (PetscBool)((some_indices[ii]%4)/2 == 1); 322 npending--; 323 if (!blockmsg) { /* Scalar stash */ 324 PetscMPIInt count; 325 if (--frame[i].pendings > 0) continue; 326 if (x->use_status) { 327 ierr = MPI_Get_count(&some_statuses[ii],intmsg ? MPIU_INT : MPIU_SCALAR,&count);CHKERRQ(ierr); 328 } else count = x->recvhdr[i].count; 329 for (j=0,recvint=frame[i].ints,recvscalar=frame[i].scalars; j<count; j++,recvint++) { 330 PetscInt loc = *recvint - X->map->rstart; 331 if (*recvint < X->map->rstart || X->map->rend <= *recvint) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Received vector entry %D out of local range [%D,%D)]",*recvint,X->map->rstart,X->map->rend); 332 switch (imode) { 333 case ADD_VALUES: 334 xarray[loc] += *recvscalar++; 335 break; 336 case INSERT_VALUES: 337 xarray[loc] = *recvscalar++; 338 break; 339 default: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Insert mode not supported 0x%x",imode); 340 } 341 } 342 } else { /* Block stash */ 343 PetscMPIInt count; 344 if (--frame[i].pendingb > 0) continue; 345 if (x->use_status) { 346 ierr = MPI_Get_count(&some_statuses[ii],intmsg ? MPIU_INT : MPIU_SCALAR,&count);CHKERRQ(ierr); 347 if (!intmsg) count /= bs; /* Convert from number of scalars to number of blocks */ 348 } else count = x->recvhdr[i].bcount; 349 for (j=0,recvint=frame[i].intb,recvscalar=frame[i].scalarb; j<count; j++,recvint++) { 350 PetscInt loc = (*recvint)*bs - X->map->rstart; 351 switch (imode) { 352 case ADD_VALUES: 353 for (k=loc; k<loc+bs; k++) xarray[k] += *recvscalar++; 354 break; 355 case INSERT_VALUES: 356 for (k=loc; k<loc+bs; k++) xarray[k] = *recvscalar++; 357 break; 358 default: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Insert mode not supported 0x%x",imode); 359 } 360 } 361 } 362 } 363 } 364 ierr = VecRestoreArray(X,&xarray);CHKERRQ(ierr); 365 ierr = MPI_Waitall(4*x->nsendranks,x->sendreqs,MPI_STATUSES_IGNORE);CHKERRQ(ierr); 366 ierr = PetscFree2(some_indices,some_statuses);CHKERRQ(ierr); 367 if (x->assembly_subset) { 368 void *dummy; /* reset segbuffers */ 369 ierr = PetscSegBufferExtractInPlace(x->segrecvint,&dummy);CHKERRQ(ierr); 370 ierr = PetscSegBufferExtractInPlace(x->segrecvscalar,&dummy);CHKERRQ(ierr); 371 } else { 372 ierr = VecAssemblyReset_MPI(X);CHKERRQ(ierr); 373 } 374 375 X->stash.insertmode = NOT_SET_VALUES; 376 X->bstash.insertmode = NOT_SET_VALUES; 377 ierr = VecStashScatterEnd_Private(&X->stash);CHKERRQ(ierr); 378 ierr = VecStashScatterEnd_Private(&X->bstash);CHKERRQ(ierr); 379 PetscFunctionReturn(0); 380 } 381 382 #undef __FUNCT__ 383 #define __FUNCT__ "VecAssemblyReset_MPI" 384 PetscErrorCode VecAssemblyReset_MPI(Vec X) 385 { 386 Vec_MPI *x = (Vec_MPI*)X->data; 387 PetscErrorCode ierr; 388 389 PetscFunctionBegin; 390 ierr = PetscFree(x->sendreqs);CHKERRQ(ierr); 391 ierr = PetscFree(x->recvreqs);CHKERRQ(ierr); 392 ierr = PetscFree(x->sendranks);CHKERRQ(ierr); 393 ierr = PetscFree(x->recvranks);CHKERRQ(ierr); 394 ierr = PetscFree(x->sendhdr);CHKERRQ(ierr); 395 ierr = PetscFree(x->recvhdr);CHKERRQ(ierr); 396 ierr = PetscFree(x->sendptrs);CHKERRQ(ierr); 397 ierr = PetscSegBufferDestroy(&x->segrecvint);CHKERRQ(ierr); 398 ierr = PetscSegBufferDestroy(&x->segrecvscalar);CHKERRQ(ierr); 399 ierr = PetscSegBufferDestroy(&x->segrecvframe);CHKERRQ(ierr); 400 PetscFunctionReturn(0); 401 } 402 403 404 #undef __FUNCT__ 405 #define __FUNCT__ "VecSetFromOptions_MPI" 406 static PetscErrorCode VecSetFromOptions_MPI(PetscOptionItems *PetscOptionsObject,Vec X) 407 { 408 PetscErrorCode ierr; 409 PetscBool flg = PETSC_FALSE,set; 410 411 PetscFunctionBegin; 412 ierr = PetscOptionsHead(PetscOptionsObject,"VecMPI Options");CHKERRQ(ierr); 413 ierr = PetscOptionsBool("-vec_assembly_bts","Use BuildTwoSided version of assembly","",flg,&flg,&set);CHKERRQ(ierr); 414 if (set) { 415 X->ops->assemblybegin = flg ? VecAssemblyBegin_MPI_BTS : VecAssemblyBegin_MPI; 416 X->ops->assemblyend = flg ? VecAssemblyEnd_MPI_BTS : VecAssemblyEnd_MPI; 417 } 418 ierr = PetscOptionsTail();CHKERRQ(ierr); 419 PetscFunctionReturn(0); 420 } 421 422 423 static struct _VecOps DvOps = { VecDuplicate_MPI, /* 1 */ 424 VecDuplicateVecs_Default, 425 VecDestroyVecs_Default, 426 VecDot_MPI, 427 VecMDot_MPI, 428 VecNorm_MPI, 429 VecTDot_MPI, 430 VecMTDot_MPI, 431 VecScale_Seq, 432 VecCopy_Seq, /* 10 */ 433 VecSet_Seq, 434 VecSwap_Seq, 435 VecAXPY_Seq, 436 VecAXPBY_Seq, 437 VecMAXPY_Seq, 438 VecAYPX_Seq, 439 VecWAXPY_Seq, 440 VecAXPBYPCZ_Seq, 441 VecPointwiseMult_Seq, 442 VecPointwiseDivide_Seq, 443 VecSetValues_MPI, /* 20 */ 444 VecAssemblyBegin_MPI, 445 VecAssemblyEnd_MPI, 446 0, 447 VecGetSize_MPI, 448 VecGetSize_Seq, 449 0, 450 VecMax_MPI, 451 VecMin_MPI, 452 VecSetRandom_Seq, 453 VecSetOption_MPI, 454 VecSetValuesBlocked_MPI, 455 VecDestroy_MPI, 456 VecView_MPI, 457 VecPlaceArray_MPI, 458 VecReplaceArray_Seq, 459 VecDot_Seq, 460 VecTDot_Seq, 461 VecNorm_Seq, 462 VecMDot_Seq, 463 VecMTDot_Seq, 464 VecLoad_Default, 465 VecReciprocal_Default, 466 VecConjugate_Seq, 467 0, 468 0, 469 VecResetArray_MPI, 470 VecSetFromOptions_MPI,/*set from options */ 471 VecMaxPointwiseDivide_Seq, 472 VecPointwiseMax_Seq, 473 VecPointwiseMaxAbs_Seq, 474 VecPointwiseMin_Seq, 475 VecGetValues_MPI, 476 0, 477 0, 478 0, 479 0, 480 0, 481 0, 482 VecStrideGather_Default, 483 VecStrideScatter_Default, 484 0, 485 0, 486 0, 487 0, 488 0, 489 VecStrideSubSetGather_Default, 490 VecStrideSubSetScatter_Default, 491 0, 492 0 493 }; 494 495 #undef __FUNCT__ 496 #define __FUNCT__ "VecCreate_MPI_Private" 497 /* 498 VecCreate_MPI_Private - Basic create routine called by VecCreate_MPI() (i.e. VecCreateMPI()), 499 VecCreateMPIWithArray(), VecCreate_Shared() (i.e. VecCreateShared()), VecCreateGhost(), 500 VecDuplicate_MPI(), VecCreateGhostWithArray(), VecDuplicate_MPI(), and VecDuplicate_Shared() 501 502 If alloc is true and array is NULL then this routine allocates the space, otherwise 503 no space is allocated. 504 */ 505 PetscErrorCode VecCreate_MPI_Private(Vec v,PetscBool alloc,PetscInt nghost,const PetscScalar array[]) 506 { 507 Vec_MPI *s; 508 PetscErrorCode ierr; 509 510 PetscFunctionBegin; 511 ierr = PetscNewLog(v,&s);CHKERRQ(ierr); 512 v->data = (void*)s; 513 ierr = PetscMemcpy(v->ops,&DvOps,sizeof(DvOps));CHKERRQ(ierr); 514 s->nghost = nghost; 515 v->petscnative = PETSC_TRUE; 516 517 ierr = PetscLayoutSetUp(v->map);CHKERRQ(ierr); 518 519 s->array = (PetscScalar*)array; 520 s->array_allocated = 0; 521 if (alloc && !array) { 522 PetscInt n = v->map->n+nghost; 523 ierr = PetscMalloc1(n,&s->array);CHKERRQ(ierr); 524 ierr = PetscLogObjectMemory((PetscObject)v,n*sizeof(PetscScalar));CHKERRQ(ierr); 525 ierr = PetscMemzero(s->array,n*sizeof(PetscScalar));CHKERRQ(ierr); 526 s->array_allocated = s->array; 527 } 528 529 /* By default parallel vectors do not have local representation */ 530 s->localrep = 0; 531 s->localupdate = 0; 532 533 v->stash.insertmode = NOT_SET_VALUES; 534 v->bstash.insertmode = NOT_SET_VALUES; 535 /* create the stashes. The block-size for bstash is set later when 536 VecSetValuesBlocked is called. 537 */ 538 ierr = VecStashCreate_Private(PetscObjectComm((PetscObject)v),1,&v->stash);CHKERRQ(ierr); 539 ierr = VecStashCreate_Private(PetscObjectComm((PetscObject)v),PetscAbs(v->map->bs),&v->bstash);CHKERRQ(ierr); 540 541 #if defined(PETSC_HAVE_MATLAB_ENGINE) 542 ierr = PetscObjectComposeFunction((PetscObject)v,"PetscMatlabEnginePut_C",VecMatlabEnginePut_Default);CHKERRQ(ierr); 543 ierr = PetscObjectComposeFunction((PetscObject)v,"PetscMatlabEngineGet_C",VecMatlabEngineGet_Default);CHKERRQ(ierr); 544 #endif 545 ierr = PetscObjectChangeTypeName((PetscObject)v,VECMPI);CHKERRQ(ierr); 546 PetscFunctionReturn(0); 547 } 548 549 /*MC 550 VECMPI - VECMPI = "mpi" - The basic parallel vector 551 552 Options Database Keys: 553 . -vec_type mpi - sets the vector type to VECMPI during a call to VecSetFromOptions() 554 555 Level: beginner 556 557 .seealso: VecCreate(), VecSetType(), VecSetFromOptions(), VecCreateMpiWithArray(), VECMPI, VecType, VecCreateMPI(), VecCreateMpi() 558 M*/ 559 560 #undef __FUNCT__ 561 #define __FUNCT__ "VecCreate_MPI" 562 PETSC_EXTERN PetscErrorCode VecCreate_MPI(Vec vv) 563 { 564 PetscErrorCode ierr; 565 566 PetscFunctionBegin; 567 ierr = VecCreate_MPI_Private(vv,PETSC_TRUE,0,0);CHKERRQ(ierr); 568 PetscFunctionReturn(0); 569 } 570 571 /*MC 572 VECSTANDARD = "standard" - A VECSEQ on one process and VECMPI on more than one process 573 574 Options Database Keys: 575 . -vec_type standard - sets a vector type to standard on calls to VecSetFromOptions() 576 577 Level: beginner 578 579 .seealso: VecCreateSeq(), VecCreateMPI() 580 M*/ 581 582 #undef __FUNCT__ 583 #define __FUNCT__ "VecCreate_Standard" 584 PETSC_EXTERN PetscErrorCode VecCreate_Standard(Vec v) 585 { 586 PetscErrorCode ierr; 587 PetscMPIInt size; 588 589 PetscFunctionBegin; 590 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)v),&size);CHKERRQ(ierr); 591 if (size == 1) { 592 ierr = VecSetType(v,VECSEQ);CHKERRQ(ierr); 593 } else { 594 ierr = VecSetType(v,VECMPI);CHKERRQ(ierr); 595 } 596 PetscFunctionReturn(0); 597 } 598 599 #undef __FUNCT__ 600 #define __FUNCT__ "VecCreateMPIWithArray" 601 /*@C 602 VecCreateMPIWithArray - Creates a parallel, array-style vector, 603 where the user provides the array space to store the vector values. 604 605 Collective on MPI_Comm 606 607 Input Parameters: 608 + comm - the MPI communicator to use 609 . bs - block size, same meaning as VecSetBlockSize() 610 . n - local vector length, cannot be PETSC_DECIDE 611 . N - global vector length (or PETSC_DECIDE to have calculated) 612 - array - the user provided array to store the vector values 613 614 Output Parameter: 615 . vv - the vector 616 617 Notes: 618 Use VecDuplicate() or VecDuplicateVecs() to form additional vectors of the 619 same type as an existing vector. 620 621 If the user-provided array is NULL, then VecPlaceArray() can be used 622 at a later stage to SET the array for storing the vector values. 623 624 PETSc does NOT free the array when the vector is destroyed via VecDestroy(). 625 The user should not free the array until the vector is destroyed. 626 627 Level: intermediate 628 629 Concepts: vectors^creating with array 630 631 .seealso: VecCreateSeqWithArray(), VecCreate(), VecDuplicate(), VecDuplicateVecs(), VecCreateGhost(), 632 VecCreateMPI(), VecCreateGhostWithArray(), VecPlaceArray() 633 634 @*/ 635 PetscErrorCode VecCreateMPIWithArray(MPI_Comm comm,PetscInt bs,PetscInt n,PetscInt N,const PetscScalar array[],Vec *vv) 636 { 637 PetscErrorCode ierr; 638 639 PetscFunctionBegin; 640 if (n == PETSC_DECIDE) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Must set local size of vector"); 641 ierr = PetscSplitOwnership(comm,&n,&N);CHKERRQ(ierr); 642 ierr = VecCreate(comm,vv);CHKERRQ(ierr); 643 ierr = VecSetSizes(*vv,n,N);CHKERRQ(ierr); 644 ierr = VecSetBlockSize(*vv,bs);CHKERRQ(ierr); 645 ierr = VecCreate_MPI_Private(*vv,PETSC_FALSE,0,array);CHKERRQ(ierr); 646 PetscFunctionReturn(0); 647 } 648 649 #undef __FUNCT__ 650 #define __FUNCT__ "VecCreateGhostWithArray" 651 /*@C 652 VecCreateGhostWithArray - Creates a parallel vector with ghost padding on each processor; 653 the caller allocates the array space. 654 655 Collective on MPI_Comm 656 657 Input Parameters: 658 + comm - the MPI communicator to use 659 . n - local vector length 660 . N - global vector length (or PETSC_DECIDE to have calculated if n is given) 661 . nghost - number of local ghost points 662 . ghosts - global indices of ghost points (or NULL if not needed), these do not need to be in increasing order (sorted) 663 - array - the space to store the vector values (as long as n + nghost) 664 665 Output Parameter: 666 . vv - the global vector representation (without ghost points as part of vector) 667 668 Notes: 669 Use VecGhostGetLocalForm() to access the local, ghosted representation 670 of the vector. 671 672 This also automatically sets the ISLocalToGlobalMapping() for this vector. 673 674 Level: advanced 675 676 Concepts: vectors^creating with array 677 Concepts: vectors^ghosted 678 679 .seealso: VecCreate(), VecGhostGetLocalForm(), VecGhostRestoreLocalForm(), 680 VecCreateGhost(), VecCreateSeqWithArray(), VecCreateMPIWithArray(), 681 VecCreateGhostBlock(), VecCreateGhostBlockWithArray(), VecMPISetGhost() 682 683 @*/ 684 PetscErrorCode VecCreateGhostWithArray(MPI_Comm comm,PetscInt n,PetscInt N,PetscInt nghost,const PetscInt ghosts[],const PetscScalar array[],Vec *vv) 685 { 686 PetscErrorCode ierr; 687 Vec_MPI *w; 688 PetscScalar *larray; 689 IS from,to; 690 ISLocalToGlobalMapping ltog; 691 PetscInt rstart,i,*indices; 692 693 PetscFunctionBegin; 694 *vv = 0; 695 696 if (n == PETSC_DECIDE) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Must set local size"); 697 if (nghost == PETSC_DECIDE) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Must set local ghost size"); 698 if (nghost < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Ghost length must be >= 0"); 699 ierr = PetscSplitOwnership(comm,&n,&N);CHKERRQ(ierr); 700 /* Create global representation */ 701 ierr = VecCreate(comm,vv);CHKERRQ(ierr); 702 ierr = VecSetSizes(*vv,n,N);CHKERRQ(ierr); 703 ierr = VecCreate_MPI_Private(*vv,PETSC_TRUE,nghost,array);CHKERRQ(ierr); 704 w = (Vec_MPI*)(*vv)->data; 705 /* Create local representation */ 706 ierr = VecGetArray(*vv,&larray);CHKERRQ(ierr); 707 ierr = VecCreateSeqWithArray(PETSC_COMM_SELF,1,n+nghost,larray,&w->localrep);CHKERRQ(ierr); 708 ierr = PetscLogObjectParent((PetscObject)*vv,(PetscObject)w->localrep);CHKERRQ(ierr); 709 ierr = VecRestoreArray(*vv,&larray);CHKERRQ(ierr); 710 711 /* 712 Create scatter context for scattering (updating) ghost values 713 */ 714 ierr = ISCreateGeneral(comm,nghost,ghosts,PETSC_COPY_VALUES,&from);CHKERRQ(ierr); 715 ierr = ISCreateStride(PETSC_COMM_SELF,nghost,n,1,&to);CHKERRQ(ierr); 716 ierr = VecScatterCreate(*vv,from,w->localrep,to,&w->localupdate);CHKERRQ(ierr); 717 ierr = PetscLogObjectParent((PetscObject)*vv,(PetscObject)w->localupdate);CHKERRQ(ierr); 718 ierr = ISDestroy(&to);CHKERRQ(ierr); 719 ierr = ISDestroy(&from);CHKERRQ(ierr); 720 721 /* set local to global mapping for ghosted vector */ 722 ierr = PetscMalloc1(n+nghost,&indices);CHKERRQ(ierr); 723 ierr = VecGetOwnershipRange(*vv,&rstart,NULL);CHKERRQ(ierr); 724 for (i=0; i<n; i++) { 725 indices[i] = rstart + i; 726 } 727 for (i=0; i<nghost; i++) { 728 indices[n+i] = ghosts[i]; 729 } 730 ierr = ISLocalToGlobalMappingCreate(comm,1,n+nghost,indices,PETSC_OWN_POINTER,<og);CHKERRQ(ierr); 731 ierr = VecSetLocalToGlobalMapping(*vv,ltog);CHKERRQ(ierr); 732 ierr = ISLocalToGlobalMappingDestroy(<og);CHKERRQ(ierr); 733 PetscFunctionReturn(0); 734 } 735 736 #undef __FUNCT__ 737 #define __FUNCT__ "VecCreateGhost" 738 /*@ 739 VecCreateGhost - Creates a parallel vector with ghost padding on each processor. 740 741 Collective on MPI_Comm 742 743 Input Parameters: 744 + comm - the MPI communicator to use 745 . n - local vector length 746 . N - global vector length (or PETSC_DECIDE to have calculated if n is given) 747 . nghost - number of local ghost points 748 - ghosts - global indices of ghost points, these do not need to be in increasing order (sorted) 749 750 Output Parameter: 751 . vv - the global vector representation (without ghost points as part of vector) 752 753 Notes: 754 Use VecGhostGetLocalForm() to access the local, ghosted representation 755 of the vector. 756 757 This also automatically sets the ISLocalToGlobalMapping() for this vector. 758 759 Level: advanced 760 761 Concepts: vectors^ghosted 762 763 .seealso: VecCreateSeq(), VecCreate(), VecDuplicate(), VecDuplicateVecs(), VecCreateMPI(), 764 VecGhostGetLocalForm(), VecGhostRestoreLocalForm(), VecGhostUpdateBegin(), 765 VecCreateGhostWithArray(), VecCreateMPIWithArray(), VecGhostUpdateEnd(), 766 VecCreateGhostBlock(), VecCreateGhostBlockWithArray(), VecMPISetGhost() 767 768 @*/ 769 PetscErrorCode VecCreateGhost(MPI_Comm comm,PetscInt n,PetscInt N,PetscInt nghost,const PetscInt ghosts[],Vec *vv) 770 { 771 PetscErrorCode ierr; 772 773 PetscFunctionBegin; 774 ierr = VecCreateGhostWithArray(comm,n,N,nghost,ghosts,0,vv);CHKERRQ(ierr); 775 PetscFunctionReturn(0); 776 } 777 778 #undef __FUNCT__ 779 #define __FUNCT__ "VecMPISetGhost" 780 /*@ 781 VecMPISetGhost - Sets the ghost points for an MPI ghost vector 782 783 Collective on Vec 784 785 Input Parameters: 786 + vv - the MPI vector 787 . nghost - number of local ghost points 788 - ghosts - global indices of ghost points, these do not need to be in increasing order (sorted) 789 790 791 Notes: 792 Use VecGhostGetLocalForm() to access the local, ghosted representation 793 of the vector. 794 795 This also automatically sets the ISLocalToGlobalMapping() for this vector. 796 797 You must call this AFTER you have set the type of the vector (with VecSetType()) and the size (with VecSetSizes()). 798 799 Level: advanced 800 801 Concepts: vectors^ghosted 802 803 .seealso: VecCreateSeq(), VecCreate(), VecDuplicate(), VecDuplicateVecs(), VecCreateMPI(), 804 VecGhostGetLocalForm(), VecGhostRestoreLocalForm(), VecGhostUpdateBegin(), 805 VecCreateGhostWithArray(), VecCreateMPIWithArray(), VecGhostUpdateEnd(), 806 VecCreateGhostBlock(), VecCreateGhostBlockWithArray() 807 808 @*/ 809 PetscErrorCode VecMPISetGhost(Vec vv,PetscInt nghost,const PetscInt ghosts[]) 810 { 811 PetscErrorCode ierr; 812 PetscBool flg; 813 814 PetscFunctionBegin; 815 ierr = PetscObjectTypeCompare((PetscObject)vv,VECMPI,&flg);CHKERRQ(ierr); 816 /* if already fully existant VECMPI then basically destroy it and rebuild with ghosting */ 817 if (flg) { 818 PetscInt n,N; 819 Vec_MPI *w; 820 PetscScalar *larray; 821 IS from,to; 822 ISLocalToGlobalMapping ltog; 823 PetscInt rstart,i,*indices; 824 MPI_Comm comm; 825 826 ierr = PetscObjectGetComm((PetscObject)vv,&comm);CHKERRQ(ierr); 827 n = vv->map->n; 828 N = vv->map->N; 829 ierr = (*vv->ops->destroy)(vv);CHKERRQ(ierr); 830 ierr = VecSetSizes(vv,n,N);CHKERRQ(ierr); 831 ierr = VecCreate_MPI_Private(vv,PETSC_TRUE,nghost,NULL);CHKERRQ(ierr); 832 w = (Vec_MPI*)(vv)->data; 833 /* Create local representation */ 834 ierr = VecGetArray(vv,&larray);CHKERRQ(ierr); 835 ierr = VecCreateSeqWithArray(PETSC_COMM_SELF,1,n+nghost,larray,&w->localrep);CHKERRQ(ierr); 836 ierr = PetscLogObjectParent((PetscObject)vv,(PetscObject)w->localrep);CHKERRQ(ierr); 837 ierr = VecRestoreArray(vv,&larray);CHKERRQ(ierr); 838 839 /* 840 Create scatter context for scattering (updating) ghost values 841 */ 842 ierr = ISCreateGeneral(comm,nghost,ghosts,PETSC_COPY_VALUES,&from);CHKERRQ(ierr); 843 ierr = ISCreateStride(PETSC_COMM_SELF,nghost,n,1,&to);CHKERRQ(ierr); 844 ierr = VecScatterCreate(vv,from,w->localrep,to,&w->localupdate);CHKERRQ(ierr); 845 ierr = PetscLogObjectParent((PetscObject)vv,(PetscObject)w->localupdate);CHKERRQ(ierr); 846 ierr = ISDestroy(&to);CHKERRQ(ierr); 847 ierr = ISDestroy(&from);CHKERRQ(ierr); 848 849 /* set local to global mapping for ghosted vector */ 850 ierr = PetscMalloc1(n+nghost,&indices);CHKERRQ(ierr); 851 ierr = VecGetOwnershipRange(vv,&rstart,NULL);CHKERRQ(ierr); 852 853 for (i=0; i<n; i++) indices[i] = rstart + i; 854 for (i=0; i<nghost; i++) indices[n+i] = ghosts[i]; 855 856 ierr = ISLocalToGlobalMappingCreate(comm,1,n+nghost,indices,PETSC_OWN_POINTER,<og);CHKERRQ(ierr); 857 ierr = VecSetLocalToGlobalMapping(vv,ltog);CHKERRQ(ierr); 858 ierr = ISLocalToGlobalMappingDestroy(<og);CHKERRQ(ierr); 859 } else if (vv->ops->create == VecCreate_MPI) SETERRQ(PetscObjectComm((PetscObject)vv),PETSC_ERR_ARG_WRONGSTATE,"Must set local or global size before setting ghosting"); 860 else if (!((PetscObject)vv)->type_name) SETERRQ(PetscObjectComm((PetscObject)vv),PETSC_ERR_ARG_WRONGSTATE,"Must set type to VECMPI before ghosting"); 861 PetscFunctionReturn(0); 862 } 863 864 865 /* ------------------------------------------------------------------------------------------*/ 866 #undef __FUNCT__ 867 #define __FUNCT__ "VecCreateGhostBlockWithArray" 868 /*@C 869 VecCreateGhostBlockWithArray - Creates a parallel vector with ghost padding on each processor; 870 the caller allocates the array space. Indices in the ghost region are based on blocks. 871 872 Collective on MPI_Comm 873 874 Input Parameters: 875 + comm - the MPI communicator to use 876 . bs - block size 877 . n - local vector length 878 . N - global vector length (or PETSC_DECIDE to have calculated if n is given) 879 . nghost - number of local ghost blocks 880 . 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) 881 - array - the space to store the vector values (as long as n + nghost*bs) 882 883 Output Parameter: 884 . vv - the global vector representation (without ghost points as part of vector) 885 886 Notes: 887 Use VecGhostGetLocalForm() to access the local, ghosted representation 888 of the vector. 889 890 n is the local vector size (total local size not the number of blocks) while nghost 891 is the number of blocks in the ghost portion, i.e. the number of elements in the ghost 892 portion is bs*nghost 893 894 Level: advanced 895 896 Concepts: vectors^creating ghosted 897 Concepts: vectors^creating with array 898 899 .seealso: VecCreate(), VecGhostGetLocalForm(), VecGhostRestoreLocalForm(), 900 VecCreateGhost(), VecCreateSeqWithArray(), VecCreateMPIWithArray(), 901 VecCreateGhostWithArray(), VecCreateGhostBlock() 902 903 @*/ 904 PetscErrorCode VecCreateGhostBlockWithArray(MPI_Comm comm,PetscInt bs,PetscInt n,PetscInt N,PetscInt nghost,const PetscInt ghosts[],const PetscScalar array[],Vec *vv) 905 { 906 PetscErrorCode ierr; 907 Vec_MPI *w; 908 PetscScalar *larray; 909 IS from,to; 910 ISLocalToGlobalMapping ltog; 911 PetscInt rstart,i,nb,*indices; 912 913 PetscFunctionBegin; 914 *vv = 0; 915 916 if (n == PETSC_DECIDE) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Must set local size"); 917 if (nghost == PETSC_DECIDE) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Must set local ghost size"); 918 if (nghost < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Ghost length must be >= 0"); 919 if (n % bs) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Local size must be a multiple of block size"); 920 ierr = PetscSplitOwnership(comm,&n,&N);CHKERRQ(ierr); 921 /* Create global representation */ 922 ierr = VecCreate(comm,vv);CHKERRQ(ierr); 923 ierr = VecSetSizes(*vv,n,N);CHKERRQ(ierr); 924 ierr = VecSetBlockSize(*vv,bs);CHKERRQ(ierr); 925 ierr = VecCreate_MPI_Private(*vv,PETSC_TRUE,nghost*bs,array);CHKERRQ(ierr); 926 w = (Vec_MPI*)(*vv)->data; 927 /* Create local representation */ 928 ierr = VecGetArray(*vv,&larray);CHKERRQ(ierr); 929 ierr = VecCreateSeqWithArray(PETSC_COMM_SELF,bs,n+bs*nghost,larray,&w->localrep);CHKERRQ(ierr); 930 ierr = PetscLogObjectParent((PetscObject)*vv,(PetscObject)w->localrep);CHKERRQ(ierr); 931 ierr = VecRestoreArray(*vv,&larray);CHKERRQ(ierr); 932 933 /* 934 Create scatter context for scattering (updating) ghost values 935 */ 936 ierr = ISCreateBlock(comm,bs,nghost,ghosts,PETSC_COPY_VALUES,&from);CHKERRQ(ierr); 937 ierr = ISCreateStride(PETSC_COMM_SELF,bs*nghost,n,1,&to);CHKERRQ(ierr); 938 ierr = VecScatterCreate(*vv,from,w->localrep,to,&w->localupdate);CHKERRQ(ierr); 939 ierr = PetscLogObjectParent((PetscObject)*vv,(PetscObject)w->localupdate);CHKERRQ(ierr); 940 ierr = ISDestroy(&to);CHKERRQ(ierr); 941 ierr = ISDestroy(&from);CHKERRQ(ierr); 942 943 /* set local to global mapping for ghosted vector */ 944 nb = n/bs; 945 ierr = PetscMalloc1(nb+nghost,&indices);CHKERRQ(ierr); 946 ierr = VecGetOwnershipRange(*vv,&rstart,NULL);CHKERRQ(ierr); 947 948 for (i=0; i<nb; i++) indices[i] = rstart + i; 949 for (i=0; i<nghost; i++) indices[nb+i] = ghosts[i]; 950 951 ierr = ISLocalToGlobalMappingCreate(comm,bs,nb+nghost,indices,PETSC_OWN_POINTER,<og);CHKERRQ(ierr); 952 ierr = VecSetLocalToGlobalMapping(*vv,ltog);CHKERRQ(ierr); 953 ierr = ISLocalToGlobalMappingDestroy(<og);CHKERRQ(ierr); 954 PetscFunctionReturn(0); 955 } 956 957 #undef __FUNCT__ 958 #define __FUNCT__ "VecCreateGhostBlock" 959 /*@ 960 VecCreateGhostBlock - Creates a parallel vector with ghost padding on each processor. 961 The indicing of the ghost points is done with blocks. 962 963 Collective on MPI_Comm 964 965 Input Parameters: 966 + comm - the MPI communicator to use 967 . bs - the block size 968 . n - local vector length 969 . N - global vector length (or PETSC_DECIDE to have calculated if n is given) 970 . nghost - number of local ghost blocks 971 - ghosts - global indices of ghost blocks, counts are by block, not by individual index, these do not need to be in increasing order (sorted) 972 973 Output Parameter: 974 . vv - the global vector representation (without ghost points as part of vector) 975 976 Notes: 977 Use VecGhostGetLocalForm() to access the local, ghosted representation 978 of the vector. 979 980 n is the local vector size (total local size not the number of blocks) while nghost 981 is the number of blocks in the ghost portion, i.e. the number of elements in the ghost 982 portion is bs*nghost 983 984 Level: advanced 985 986 Concepts: vectors^ghosted 987 988 .seealso: VecCreateSeq(), VecCreate(), VecDuplicate(), VecDuplicateVecs(), VecCreateMPI(), 989 VecGhostGetLocalForm(), VecGhostRestoreLocalForm(), 990 VecCreateGhostWithArray(), VecCreateMPIWithArray(), VecCreateGhostBlockWithArray() 991 992 @*/ 993 PetscErrorCode VecCreateGhostBlock(MPI_Comm comm,PetscInt bs,PetscInt n,PetscInt N,PetscInt nghost,const PetscInt ghosts[],Vec *vv) 994 { 995 PetscErrorCode ierr; 996 997 PetscFunctionBegin; 998 ierr = VecCreateGhostBlockWithArray(comm,bs,n,N,nghost,ghosts,0,vv);CHKERRQ(ierr); 999 PetscFunctionReturn(0); 1000 } 1001