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