Actual source code: dmredundant.c

petsc-3.5.1 2014-08-06
Report Typos and Errors
  1: #include <petsc-private/dmimpl.h>
  2: #include <petscdmredundant.h>   /*I      "petscdmredundant.h" I*/

  4: typedef struct  {
  5:   PetscInt rank;                /* owner */
  6:   PetscInt N;                   /* total number of dofs */
  7:   PetscInt n;                   /* owned number of dofs, n=N on owner, n=0 on non-owners */
  8: } DM_Redundant;

 12: static PetscErrorCode DMCreateMatrix_Redundant(DM dm,Mat *J)
 13: {
 14:   DM_Redundant           *red = (DM_Redundant*)dm->data;
 15:   PetscErrorCode         ierr;
 16:   ISLocalToGlobalMapping ltog;
 17:   PetscInt               i,rstart,rend,*cols;
 18:   PetscScalar            *vals;

 21:   MatCreate(PetscObjectComm((PetscObject)dm),J);
 22:   MatSetSizes(*J,red->n,red->n,red->N,red->N);
 23:   MatSetType(*J,dm->mattype);
 24:   MatSeqAIJSetPreallocation(*J,red->n,NULL);
 25:   MatSeqBAIJSetPreallocation(*J,1,red->n,NULL);
 26:   MatMPIAIJSetPreallocation(*J,red->n,NULL,red->N-red->n,NULL);
 27:   MatMPIBAIJSetPreallocation(*J,1,red->n,NULL,red->N-red->n,NULL);

 29:   DMGetLocalToGlobalMapping(dm,&ltog);
 30:   MatSetLocalToGlobalMapping(*J,ltog,ltog);

 32:   PetscMalloc2(red->N,&cols,red->N,&vals);
 33:   for (i=0; i<red->N; i++) {
 34:     cols[i] = i;
 35:     vals[i] = 0.0;
 36:   }
 37:   MatGetOwnershipRange(*J,&rstart,&rend);
 38:   for (i=rstart; i<rend; i++) {
 39:     MatSetValues(*J,1,&i,red->N,cols,vals,INSERT_VALUES);
 40:   }
 41:   PetscFree2(cols,vals);
 42:   MatAssemblyBegin(*J,MAT_FINAL_ASSEMBLY);
 43:   MatAssemblyEnd(*J,MAT_FINAL_ASSEMBLY);
 44:   return(0);
 45: }

 49: static PetscErrorCode DMDestroy_Redundant(DM dm)
 50: {

 54:   PetscObjectComposeFunction((PetscObject)dm,"DMRedundantSetSize_C",NULL);
 55:   PetscObjectComposeFunction((PetscObject)dm,"DMRedundantGetSize_C",NULL);
 56:   /* This was originally freed in DMDestroy(), but that prevents reference counting of backend objects */
 57:   PetscFree(dm->data);
 58:   return(0);
 59: }

 63: static PetscErrorCode DMCreateGlobalVector_Redundant(DM dm,Vec *gvec)
 64: {
 65:   PetscErrorCode         ierr;
 66:   DM_Redundant           *red = (DM_Redundant*)dm->data;
 67:   ISLocalToGlobalMapping ltog;

 72:   *gvec = 0;
 73:   VecCreate(PetscObjectComm((PetscObject)dm),gvec);
 74:   VecSetSizes(*gvec,red->n,red->N);
 75:   VecSetType(*gvec,dm->vectype);
 76:   DMGetLocalToGlobalMapping(dm,&ltog);
 77:   VecSetLocalToGlobalMapping(*gvec,ltog);
 78:   VecSetDM(*gvec,dm);
 79:   return(0);
 80: }

 84: static PetscErrorCode DMCreateLocalVector_Redundant(DM dm,Vec *lvec)
 85: {
 87:   DM_Redundant   *red = (DM_Redundant*)dm->data;

 92:   *lvec = 0;
 93:   VecCreate(PETSC_COMM_SELF,lvec);
 94:   VecSetSizes(*lvec,red->N,red->N);
 95:   VecSetType(*lvec,dm->vectype);
 96:   VecSetDM(*lvec,dm);
 97:   return(0);
 98: }

102: static PetscErrorCode DMLocalToGlobalBegin_Redundant(DM dm,Vec l,InsertMode imode,Vec g)
103: {
104:   PetscErrorCode    ierr;
105:   DM_Redundant      *red = (DM_Redundant*)dm->data;
106:   const PetscScalar *lv;
107:   PetscScalar       *gv;
108:   PetscMPIInt       rank;

111:   MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&rank);
112:   VecGetArrayRead(l,&lv);
113:   VecGetArray(g,&gv);
114:   switch (imode) {
115:   case ADD_VALUES:
116:   case MAX_VALUES:
117:   {
118:     void        *source;
119:     PetscScalar *buffer;
120:     PetscInt    i;
121:     if (rank == red->rank) {
122: #if defined(PETSC_HAVE_MPI_IN_PLACE)
123:       buffer = gv;
124:       source = MPI_IN_PLACE;
125: #else
126:       PetscMalloc1(red->N,&buffer);
127:       source = buffer;
128: #endif
129:       if (imode == ADD_VALUES) for (i=0; i<red->N; i++) buffer[i] = gv[i] + lv[i];
130: #if !defined(PETSC_USE_COMPLEX)
131:       if (imode == MAX_VALUES) for (i=0; i<red->N; i++) buffer[i] = PetscMax(gv[i],lv[i]);
132: #endif
133:     } else source = (void*)lv;
134:     MPI_Reduce(source,gv,red->N,MPIU_SCALAR,(imode == ADD_VALUES) ? MPIU_SUM : MPIU_MAX,red->rank,PetscObjectComm((PetscObject)dm));
135: #if !defined(PETSC_HAVE_MPI_IN_PLACE)
136:     if (rank == red->rank) {PetscFree(buffer);}
137: #endif
138:   } break;
139:   case INSERT_VALUES:
140:     PetscMemcpy(gv,lv,red->n*sizeof(PetscScalar));
141:     break;
142:   default: SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"InsertMode not supported");
143:   }
144:   VecRestoreArrayRead(l,&lv);
145:   VecRestoreArray(g,&gv);
146:   return(0);
147: }

151: static PetscErrorCode DMLocalToGlobalEnd_Redundant(DM dm,Vec l,InsertMode imode,Vec g)
152: {
154:   return(0);
155: }

159: static PetscErrorCode DMGlobalToLocalBegin_Redundant(DM dm,Vec g,InsertMode imode,Vec l)
160: {
161:   PetscErrorCode    ierr;
162:   DM_Redundant      *red = (DM_Redundant*)dm->data;
163:   const PetscScalar *gv;
164:   PetscScalar       *lv;

167:   VecGetArrayRead(g,&gv);
168:   VecGetArray(l,&lv);
169:   switch (imode) {
170:   case INSERT_VALUES:
171:     if (red->n) {PetscMemcpy(lv,gv,red->n*sizeof(PetscScalar));}
172:     MPI_Bcast(lv,red->N,MPIU_SCALAR,red->rank,PetscObjectComm((PetscObject)dm));
173:     break;
174:   default: SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"InsertMode not supported");
175:   }
176:   VecRestoreArrayRead(g,&gv);
177:   VecRestoreArray(l,&lv);
178:   return(0);
179: }

183: static PetscErrorCode DMGlobalToLocalEnd_Redundant(DM dm,Vec g,InsertMode imode,Vec l)
184: {
186:   return(0);
187: }

191: static PetscErrorCode DMSetUp_Redundant(DM dm)
192: {
194:   DM_Redundant   *red = (DM_Redundant*)dm->data;
195:   PetscInt       i,*globals;

198:   PetscMalloc1(red->N,&globals);
199:   for (i=0; i<red->N; i++) globals[i] = i;
200:   ISLocalToGlobalMappingCreate(PETSC_COMM_SELF,1,red->N,globals,PETSC_OWN_POINTER,&dm->ltogmap);
201:   return(0);
202: }

206: static PetscErrorCode DMView_Redundant(DM dm,PetscViewer viewer)
207: {
209:   DM_Redundant   *red = (DM_Redundant*)dm->data;
210:   PetscBool      iascii;

213:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
214:   if (iascii) {
215:     PetscViewerASCIIPrintf(viewer,"redundant: rank=%D N=%D\n",red->rank,red->N);
216:   }
217:   return(0);
218: }

222: static PetscErrorCode DMCreateColoring_Redundant(DM dm,ISColoringType ctype,ISColoring *coloring)
223: {
224:   DM_Redundant    *red = (DM_Redundant*)dm->data;
225:   PetscErrorCode  ierr;
226:   PetscInt        i,nloc;
227:   ISColoringValue *colors;

230:   switch (ctype) {
231:   case IS_COLORING_GLOBAL:
232:     nloc = red->n;
233:     break;
234:   case IS_COLORING_GHOSTED:
235:     nloc = red->N;
236:     break;
237:   default: SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONG,"Unknown ISColoringType %d",(int)ctype);
238:   }
239:   PetscMalloc1(nloc,&colors);
240:   for (i=0; i<nloc; i++) colors[i] = i;
241:   ISColoringCreate(PetscObjectComm((PetscObject)dm),red->N,nloc,colors,coloring);
242:   ISColoringSetType(*coloring,ctype);
243:   return(0);
244: }

248: static PetscErrorCode DMRefine_Redundant(DM dmc,MPI_Comm comm,DM *dmf)
249: {
251:   PetscMPIInt    flag;
252:   DM_Redundant   *redc = (DM_Redundant*)dmc->data;

255:   if (comm == MPI_COMM_NULL) {
256:     PetscObjectGetComm((PetscObject)dmc,&comm);
257:   }
258:   MPI_Comm_compare(PetscObjectComm((PetscObject)dmc),comm,&flag);
259:   if (flag != MPI_CONGRUENT && flag != MPI_IDENT) SETERRQ(PetscObjectComm((PetscObject)dmc),PETSC_ERR_SUP,"cannot change communicators");
260:   DMRedundantCreate(comm,redc->rank,redc->N,dmf);
261:   return(0);
262: }

266: static PetscErrorCode DMCoarsen_Redundant(DM dmf,MPI_Comm comm,DM *dmc)
267: {
269:   PetscMPIInt    flag;
270:   DM_Redundant   *redf = (DM_Redundant*)dmf->data;

273:   if (comm == MPI_COMM_NULL) {
274:     PetscObjectGetComm((PetscObject)dmf,&comm);
275:   }
276:   MPI_Comm_compare(PetscObjectComm((PetscObject)dmf),comm,&flag);
277:   if (flag != MPI_CONGRUENT && flag != MPI_IDENT) SETERRQ(PetscObjectComm((PetscObject)dmf),PETSC_ERR_SUP,"cannot change communicators");
278:   DMRedundantCreate(comm,redf->rank,redf->N,dmc);
279:   return(0);
280: }

284: static PetscErrorCode DMCreateInterpolation_Redundant(DM dmc,DM dmf,Mat *P,Vec *scale)
285: {
287:   DM_Redundant   *redc = (DM_Redundant*)dmc->data;
288:   DM_Redundant   *redf = (DM_Redundant*)dmf->data;
289:   PetscMPIInt    flag;
290:   PetscInt       i,rstart,rend;

293:   MPI_Comm_compare(PetscObjectComm((PetscObject)dmc),PetscObjectComm((PetscObject)dmf),&flag);
294:   if (flag != MPI_CONGRUENT && flag != MPI_IDENT) SETERRQ(PetscObjectComm((PetscObject)dmf),PETSC_ERR_SUP,"cannot change communicators");
295:   if (redc->rank != redf->rank) SETERRQ(PetscObjectComm((PetscObject)dmf),PETSC_ERR_ARG_INCOMP,"Owning rank does not match");
296:   if (redc->N != redf->N) SETERRQ(PetscObjectComm((PetscObject)dmf),PETSC_ERR_ARG_INCOMP,"Global size does not match");
297:   MatCreate(PetscObjectComm((PetscObject)dmc),P);
298:   MatSetSizes(*P,redc->n,redc->n,redc->N,redc->N);
299:   MatSetType(*P,MATAIJ);
300:   MatSeqAIJSetPreallocation(*P,1,0);
301:   MatMPIAIJSetPreallocation(*P,1,0,0,0);
302:   MatGetOwnershipRange(*P,&rstart,&rend);
303:   for (i=rstart; i<rend; i++) {MatSetValue(*P,i,i,1.0,INSERT_VALUES);}
304:   MatAssemblyBegin(*P,MAT_FINAL_ASSEMBLY);
305:   MatAssemblyEnd(*P,MAT_FINAL_ASSEMBLY);
306:   if (scale) {DMCreateInterpolationScale(dmc,dmf,*P,scale);}
307:   return(0);
308: }

312: /*@
313:     DMRedundantSetSize - Sets the size of a densely coupled redundant object

315:     Collective on DM

317:     Input Parameter:
318: +   dm - redundant DM
319: .   rank - rank of process to own redundant degrees of freedom
320: -   N - total number of redundant degrees of freedom

322:     Level: advanced

324: .seealso DMDestroy(), DMCreateGlobalVector(), DMRedundantCreate(), DMRedundantGetSize()
325: @*/
326: PetscErrorCode DMRedundantSetSize(DM dm,PetscInt rank,PetscInt N)
327: {

335:   PetscTryMethod(dm,"DMRedundantSetSize_C",(DM,PetscInt,PetscInt),(dm,rank,N));
336:   return(0);
337: }

341: /*@
342:     DMRedundantGetSize - Gets the size of a densely coupled redundant object

344:     Not Collective

346:     Input Parameter:
347: +   dm - redundant DM

349:     Output Parameters:
350: +   rank - rank of process to own redundant degrees of freedom (or NULL)
351: -   N - total number of redundant degrees of freedom (or NULL)

353:     Level: advanced

355: .seealso DMDestroy(), DMCreateGlobalVector(), DMRedundantCreate(), DMRedundantSetSize()
356: @*/
357: PetscErrorCode DMRedundantGetSize(DM dm,PetscInt *rank,PetscInt *N)
358: {

364:   PetscUseMethod(dm,"DMRedundantGetSize_C",(DM,PetscInt*,PetscInt*),(dm,rank,N));
365:   return(0);
366: }

370: static PetscErrorCode DMRedundantSetSize_Redundant(DM dm,PetscInt rank,PetscInt N)
371: {
372:   DM_Redundant   *red = (DM_Redundant*)dm->data;
374:   PetscMPIInt    myrank;

377:   MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&myrank);
378:   red->rank = rank;
379:   red->N    = N;
380:   red->n    = (myrank == rank) ? N : 0;
381:   return(0);
382: }

386: static PetscErrorCode DMRedundantGetSize_Redundant(DM dm,PetscInt *rank,PetscInt *N)
387: {
388:   DM_Redundant *red = (DM_Redundant*)dm->data;

391:   if (rank) *rank = red->rank;
392:   if (N)    *N = red->N;
393:   return(0);
394: }

396: /*MC
397:    DMREDUNDANT = "redundant" - A DM object that is used to manage data for a small set of dense globally coupled variables.
398:          In the global representation of the vector the variables are all stored on a single MPI process (all the other MPI processes have
399:          no variables) in the local representation all the variables are stored on ALL the MPI processes (because they are all needed for each
400:          processes local computations).

402:          This DM is generally used inside a DMCOMPOSITE object. For example, it may be used to store continuation parameters for a bifurcation problem.


405:   Level: intermediate

407: .seealso: DMType, DMCOMPOSITE, DMCreateRedundant(), DMCreate(), DMRedundantSetSize(), DMRedundantGetSize()
408: M*/

412: PETSC_EXTERN PetscErrorCode DMCreate_Redundant(DM dm)
413: {
415:   DM_Redundant   *red;

418:   PetscNewLog(dm,&red);
419:   dm->data = red;

421:   PetscObjectChangeTypeName((PetscObject)dm,DMREDUNDANT);

423:   dm->ops->setup              = DMSetUp_Redundant;
424:   dm->ops->view               = DMView_Redundant;
425:   dm->ops->createglobalvector = DMCreateGlobalVector_Redundant;
426:   dm->ops->createlocalvector  = DMCreateLocalVector_Redundant;
427:   dm->ops->creatematrix       = DMCreateMatrix_Redundant;
428:   dm->ops->destroy            = DMDestroy_Redundant;
429:   dm->ops->globaltolocalbegin = DMGlobalToLocalBegin_Redundant;
430:   dm->ops->globaltolocalend   = DMGlobalToLocalEnd_Redundant;
431:   dm->ops->localtoglobalbegin = DMLocalToGlobalBegin_Redundant;
432:   dm->ops->localtoglobalend   = DMLocalToGlobalEnd_Redundant;
433:   dm->ops->refine             = DMRefine_Redundant;
434:   dm->ops->coarsen            = DMCoarsen_Redundant;
435:   dm->ops->createinterpolation= DMCreateInterpolation_Redundant;
436:   dm->ops->getcoloring        = DMCreateColoring_Redundant;

438:   PetscObjectComposeFunction((PetscObject)dm,"DMRedundantSetSize_C",DMRedundantSetSize_Redundant);
439:   PetscObjectComposeFunction((PetscObject)dm,"DMRedundantGetSize_C",DMRedundantGetSize_Redundant);
440:   return(0);
441: }

445: /*@C
446:     DMRedundantCreate - Creates a DM object, used to manage data for dense globally coupled variables

448:     Collective on MPI_Comm

450:     Input Parameter:
451: +   comm - the processors that will share the global vector
452: .   rank - rank to own the redundant values
453: -   N - total number of degrees of freedom

455:     Output Parameters:
456: .   red - the redundant DM

458:     Level: advanced

460: .seealso DMDestroy(), DMCreateGlobalVector(), DMCreateMatrix(), DMCompositeAddDM(), DMREDUNDANT, DMSetType(), DMRedundantSetSize(), DMRedundantGetSize()

462: @*/
463: PetscErrorCode DMRedundantCreate(MPI_Comm comm,PetscInt rank,PetscInt N,DM *dm)
464: {

469:   DMCreate(comm,dm);
470:   DMSetType(*dm,DMREDUNDANT);
471:   DMRedundantSetSize(*dm,rank,N);
472:   DMSetUp(*dm);
473:   return(0);
474: }