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