107475bc1SBarry Smith #include <petsc-private/dmimpl.h> 28ac4e037SJed Brown #include <petscdmredundant.h> /*I "petscdmredundant.h" I*/ 38ac4e037SJed Brown 48ac4e037SJed Brown typedef struct { 58ac4e037SJed Brown PetscInt rank; /* owner */ 68ac4e037SJed Brown PetscInt N; /* total number of dofs */ 78ac4e037SJed Brown PetscInt n; /* owned number of dofs, n=N on owner, n=0 on non-owners */ 88ac4e037SJed Brown } DM_Redundant; 98ac4e037SJed Brown 108ac4e037SJed Brown #undef __FUNCT__ 11950540a4SJed Brown #define __FUNCT__ "DMCreateMatrix_Redundant" 1219fd82e9SBarry Smith static PetscErrorCode DMCreateMatrix_Redundant(DM dm,MatType mtype,Mat *J) 138ac4e037SJed Brown { 148ac4e037SJed Brown DM_Redundant *red = (DM_Redundant*)dm->data; 158ac4e037SJed Brown PetscErrorCode ierr; 168ac4e037SJed Brown ISLocalToGlobalMapping ltog,ltogb; 17d2e2b695SJed Brown PetscInt i,rstart,rend,*cols; 18d2e2b695SJed Brown PetscScalar *vals; 198ac4e037SJed Brown 208ac4e037SJed Brown PetscFunctionBegin; 21ce94432eSBarry Smith ierr = MatCreate(PetscObjectComm((PetscObject)dm),J);CHKERRQ(ierr); 228ac4e037SJed Brown ierr = MatSetSizes(*J,red->n,red->n,red->N,red->N);CHKERRQ(ierr); 238ac4e037SJed Brown ierr = MatSetType(*J,mtype);CHKERRQ(ierr); 240298fd71SBarry Smith ierr = MatSeqAIJSetPreallocation(*J,red->n,NULL);CHKERRQ(ierr); 250298fd71SBarry Smith ierr = MatSeqBAIJSetPreallocation(*J,1,red->n,NULL);CHKERRQ(ierr); 260298fd71SBarry Smith ierr = MatMPIAIJSetPreallocation(*J,red->n,NULL,red->N-red->n,NULL);CHKERRQ(ierr); 270298fd71SBarry Smith ierr = MatMPIBAIJSetPreallocation(*J,1,red->n,NULL,red->N-red->n,NULL);CHKERRQ(ierr); 288ac4e037SJed Brown 298ac4e037SJed Brown ierr = DMGetLocalToGlobalMapping(dm,<og);CHKERRQ(ierr); 308ac4e037SJed Brown ierr = DMGetLocalToGlobalMappingBlock(dm,<ogb);CHKERRQ(ierr); 318ac4e037SJed Brown ierr = MatSetLocalToGlobalMapping(*J,ltog,ltog);CHKERRQ(ierr); 328ac4e037SJed Brown ierr = MatSetLocalToGlobalMappingBlock(*J,ltogb,ltogb);CHKERRQ(ierr); 33d2e2b695SJed Brown 34d2e2b695SJed Brown ierr = PetscMalloc2(red->N,PetscInt,&cols,red->N,PetscScalar,&vals);CHKERRQ(ierr); 35d2e2b695SJed Brown for (i=0; i<red->N; i++) { 36d2e2b695SJed Brown cols[i] = i; 37d2e2b695SJed Brown vals[i] = 0.0; 38d2e2b695SJed Brown } 39d2e2b695SJed Brown ierr = MatGetOwnershipRange(*J,&rstart,&rend);CHKERRQ(ierr); 40d2e2b695SJed Brown for (i=rstart; i<rend; i++) { 41d2e2b695SJed Brown ierr = MatSetValues(*J,1,&i,red->N,cols,vals,INSERT_VALUES);CHKERRQ(ierr); 42d2e2b695SJed Brown } 43d2e2b695SJed Brown ierr = PetscFree2(cols,vals);CHKERRQ(ierr); 44d2e2b695SJed Brown ierr = MatAssemblyBegin(*J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 45d2e2b695SJed Brown ierr = MatAssemblyEnd(*J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 468ac4e037SJed Brown PetscFunctionReturn(0); 478ac4e037SJed Brown } 488ac4e037SJed Brown 498ac4e037SJed Brown #undef __FUNCT__ 508ac4e037SJed Brown #define __FUNCT__ "DMDestroy_Redundant" 518ac4e037SJed Brown static PetscErrorCode DMDestroy_Redundant(DM dm) 528ac4e037SJed Brown { 538ac4e037SJed Brown PetscErrorCode ierr; 548ac4e037SJed Brown 558ac4e037SJed Brown PetscFunctionBegin; 5600de8ff0SBarry Smith ierr = PetscObjectComposeFunction((PetscObject)dm,"DMRedundantSetSize_C","",NULL);CHKERRQ(ierr); 5700de8ff0SBarry Smith ierr = PetscObjectComposeFunction((PetscObject)dm,"DMRedundantGetSize_C","",NULL);CHKERRQ(ierr); 58435a35e8SMatthew G Knepley /* This was originally freed in DMDestroy(), but that prevents reference counting of backend objects */ 59435a35e8SMatthew G Knepley ierr = PetscFree(dm->data);CHKERRQ(ierr); 608ac4e037SJed Brown PetscFunctionReturn(0); 618ac4e037SJed Brown } 628ac4e037SJed Brown 638ac4e037SJed Brown #undef __FUNCT__ 648ac4e037SJed Brown #define __FUNCT__ "DMCreateGlobalVector_Redundant" 658ac4e037SJed Brown static PetscErrorCode DMCreateGlobalVector_Redundant(DM dm,Vec *gvec) 668ac4e037SJed Brown { 678ac4e037SJed Brown PetscErrorCode ierr; 688ac4e037SJed Brown DM_Redundant *red = (DM_Redundant*)dm->data; 698ac4e037SJed Brown ISLocalToGlobalMapping ltog,ltogb; 708ac4e037SJed Brown 718ac4e037SJed Brown PetscFunctionBegin; 728ac4e037SJed Brown PetscValidHeaderSpecific(dm,DM_CLASSID,1); 738ac4e037SJed Brown PetscValidPointer(gvec,2); 748ac4e037SJed Brown *gvec = 0; 75ce94432eSBarry Smith ierr = VecCreate(PetscObjectComm((PetscObject)dm),gvec);CHKERRQ(ierr); 768ac4e037SJed Brown ierr = VecSetSizes(*gvec,red->n,red->N);CHKERRQ(ierr); 778ac4e037SJed Brown ierr = VecSetType(*gvec,dm->vectype);CHKERRQ(ierr); 788ac4e037SJed Brown ierr = DMGetLocalToGlobalMapping(dm,<og);CHKERRQ(ierr); 798ac4e037SJed Brown ierr = DMGetLocalToGlobalMappingBlock(dm,<ogb);CHKERRQ(ierr); 808ac4e037SJed Brown ierr = VecSetLocalToGlobalMapping(*gvec,ltog);CHKERRQ(ierr); 818ac4e037SJed Brown ierr = VecSetLocalToGlobalMappingBlock(*gvec,ltog);CHKERRQ(ierr); 82c688c046SMatthew G Knepley ierr = VecSetDM(*gvec,dm);CHKERRQ(ierr); 838ac4e037SJed Brown PetscFunctionReturn(0); 848ac4e037SJed Brown } 858ac4e037SJed Brown 868ac4e037SJed Brown #undef __FUNCT__ 878ac4e037SJed Brown #define __FUNCT__ "DMCreateLocalVector_Redundant" 888ac4e037SJed Brown static PetscErrorCode DMCreateLocalVector_Redundant(DM dm,Vec *lvec) 898ac4e037SJed Brown { 908ac4e037SJed Brown PetscErrorCode ierr; 918ac4e037SJed Brown DM_Redundant *red = (DM_Redundant*)dm->data; 928ac4e037SJed Brown 938ac4e037SJed Brown PetscFunctionBegin; 948ac4e037SJed Brown PetscValidHeaderSpecific(dm,DM_CLASSID,1); 958ac4e037SJed Brown PetscValidPointer(lvec,2); 968ac4e037SJed Brown *lvec = 0; 978ac4e037SJed Brown ierr = VecCreate(PETSC_COMM_SELF,lvec);CHKERRQ(ierr); 988ac4e037SJed Brown ierr = VecSetSizes(*lvec,red->N,red->N);CHKERRQ(ierr); 998ac4e037SJed Brown ierr = VecSetType(*lvec,dm->vectype);CHKERRQ(ierr); 100c688c046SMatthew G Knepley ierr = VecSetDM(*lvec,dm);CHKERRQ(ierr); 1018ac4e037SJed Brown PetscFunctionReturn(0); 1028ac4e037SJed Brown } 1038ac4e037SJed Brown 1048ac4e037SJed Brown #undef __FUNCT__ 1058ac4e037SJed Brown #define __FUNCT__ "DMLocalToGlobalBegin_Redundant" 1068ac4e037SJed Brown static PetscErrorCode DMLocalToGlobalBegin_Redundant(DM dm,Vec l,InsertMode imode,Vec g) 1078ac4e037SJed Brown { 1088ac4e037SJed Brown PetscErrorCode ierr; 1098ac4e037SJed Brown DM_Redundant *red = (DM_Redundant*)dm->data; 1108ac4e037SJed Brown const PetscScalar *lv; 1118ac4e037SJed Brown PetscScalar *gv; 1128ac4e037SJed Brown PetscMPIInt rank; 1138ac4e037SJed Brown 1148ac4e037SJed Brown PetscFunctionBegin; 115ce94432eSBarry Smith ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&rank);CHKERRQ(ierr); 1168ac4e037SJed Brown ierr = VecGetArrayRead(l,&lv);CHKERRQ(ierr); 1178ac4e037SJed Brown ierr = VecGetArray(g,&gv);CHKERRQ(ierr); 1188ac4e037SJed Brown switch (imode) { 1198ac4e037SJed Brown case ADD_VALUES: 1208ac4e037SJed Brown case MAX_VALUES: 1218ac4e037SJed Brown { 1228ac4e037SJed Brown void *source; 1238ac4e037SJed Brown PetscScalar *buffer; 1248ac4e037SJed Brown PetscInt i; 1258ac4e037SJed Brown if (rank == red->rank) { 1268ac4e037SJed Brown #if defined(PETSC_HAVE_MPI_IN_PLACE) 1278ac4e037SJed Brown buffer = gv; 1288ac4e037SJed Brown source = MPI_IN_PLACE; 1298ac4e037SJed Brown #else 1308ac4e037SJed Brown ierr = PetscMalloc(red->N*sizeof(PetscScalar),&buffer);CHKERRQ(ierr); 1318ac4e037SJed Brown source = buffer; 1328ac4e037SJed Brown #endif 1338ac4e037SJed Brown if (imode == ADD_VALUES) for (i=0; i<red->N; i++) buffer[i] = gv[i] + lv[i]; 1348ac4e037SJed Brown #if !defined(PETSC_USE_COMPLEX) 1358ac4e037SJed Brown if (imode == MAX_VALUES) for (i=0; i<red->N; i++) buffer[i] = PetscMax(gv[i],lv[i]); 1368ac4e037SJed Brown #endif 1378865f1eaSKarl Rupp } else source = (void*)lv; 138ce94432eSBarry Smith ierr = MPI_Reduce(source,gv,red->N,MPIU_SCALAR,(imode == ADD_VALUES) ? MPI_SUM : MPI_MAX,red->rank,PetscObjectComm((PetscObject)dm));CHKERRQ(ierr); 1398ac4e037SJed Brown #if !defined(PETSC_HAVE_MPI_IN_PLACE) 14099839c1eSBarry Smith if (rank == red->rank) {ierr = PetscFree(buffer);CHKERRQ(ierr);} 1418ac4e037SJed Brown #endif 1428ac4e037SJed Brown } break; 1438ac4e037SJed Brown case INSERT_VALUES: 1448ac4e037SJed Brown ierr = PetscMemcpy(gv,lv,red->n*sizeof(PetscScalar));CHKERRQ(ierr); 1458ac4e037SJed Brown break; 146ce94432eSBarry Smith default: SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"InsertMode not supported"); 1478ac4e037SJed Brown } 1488ac4e037SJed Brown ierr = VecRestoreArrayRead(l,&lv);CHKERRQ(ierr); 1498ac4e037SJed Brown ierr = VecRestoreArray(g,&gv);CHKERRQ(ierr); 1508ac4e037SJed Brown PetscFunctionReturn(0); 1518ac4e037SJed Brown } 1528ac4e037SJed Brown 1538ac4e037SJed Brown #undef __FUNCT__ 1548ac4e037SJed Brown #define __FUNCT__ "DMLocalToGlobalEnd_Redundant" 1558ac4e037SJed Brown static PetscErrorCode DMLocalToGlobalEnd_Redundant(DM dm,Vec l,InsertMode imode,Vec g) 1568ac4e037SJed Brown { 1578ac4e037SJed Brown PetscFunctionBegin; 1588ac4e037SJed Brown PetscFunctionReturn(0); 1598ac4e037SJed Brown } 1608ac4e037SJed Brown 1618ac4e037SJed Brown #undef __FUNCT__ 1628ac4e037SJed Brown #define __FUNCT__ "DMGlobalToLocalBegin_Redundant" 1638ac4e037SJed Brown static PetscErrorCode DMGlobalToLocalBegin_Redundant(DM dm,Vec g,InsertMode imode,Vec l) 1648ac4e037SJed Brown { 1658ac4e037SJed Brown PetscErrorCode ierr; 1668ac4e037SJed Brown DM_Redundant *red = (DM_Redundant*)dm->data; 1678ac4e037SJed Brown const PetscScalar *gv; 1688ac4e037SJed Brown PetscScalar *lv; 1698ac4e037SJed Brown 1708ac4e037SJed Brown PetscFunctionBegin; 1718ac4e037SJed Brown ierr = VecGetArrayRead(g,&gv);CHKERRQ(ierr); 1728ac4e037SJed Brown ierr = VecGetArray(l,&lv);CHKERRQ(ierr); 1738ac4e037SJed Brown switch (imode) { 1748ac4e037SJed Brown case INSERT_VALUES: 1758ac4e037SJed Brown if (red->n) {ierr = PetscMemcpy(lv,gv,red->n*sizeof(PetscScalar));CHKERRQ(ierr);} 176ce94432eSBarry Smith ierr = MPI_Bcast(lv,red->N,MPIU_SCALAR,red->rank,PetscObjectComm((PetscObject)dm));CHKERRQ(ierr); 1778ac4e037SJed Brown break; 178ce94432eSBarry Smith default: SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"InsertMode not supported"); 1798ac4e037SJed Brown } 1808ac4e037SJed Brown ierr = VecRestoreArrayRead(g,&gv);CHKERRQ(ierr); 1818ac4e037SJed Brown ierr = VecRestoreArray(l,&lv);CHKERRQ(ierr); 1828ac4e037SJed Brown PetscFunctionReturn(0); 1838ac4e037SJed Brown } 1848ac4e037SJed Brown 1858ac4e037SJed Brown #undef __FUNCT__ 1868ac4e037SJed Brown #define __FUNCT__ "DMGlobalToLocalEnd_Redundant" 1878ac4e037SJed Brown static PetscErrorCode DMGlobalToLocalEnd_Redundant(DM dm,Vec g,InsertMode imode,Vec l) 1888ac4e037SJed Brown { 1898ac4e037SJed Brown PetscFunctionBegin; 1908ac4e037SJed Brown PetscFunctionReturn(0); 1918ac4e037SJed Brown } 1928ac4e037SJed Brown 1938ac4e037SJed Brown #undef __FUNCT__ 1948ac4e037SJed Brown #define __FUNCT__ "DMSetUp_Redundant" 1958ac4e037SJed Brown static PetscErrorCode DMSetUp_Redundant(DM dm) 1968ac4e037SJed Brown { 1978ac4e037SJed Brown PetscErrorCode ierr; 1988ac4e037SJed Brown DM_Redundant *red = (DM_Redundant*)dm->data; 1998ac4e037SJed Brown PetscInt i,*globals; 2008ac4e037SJed Brown 2018ac4e037SJed Brown PetscFunctionBegin; 2028ac4e037SJed Brown ierr = PetscMalloc(red->N*sizeof(PetscInt),&globals);CHKERRQ(ierr); 2038ac4e037SJed Brown for (i=0; i<red->N; i++) globals[i] = i; 2048ac4e037SJed Brown ierr = ISLocalToGlobalMappingCreate(PETSC_COMM_SELF,red->N,globals,PETSC_OWN_POINTER,&dm->ltogmap);CHKERRQ(ierr); 2058ac4e037SJed Brown ierr = PetscObjectReference((PetscObject)dm->ltogmap);CHKERRQ(ierr); 2068ac4e037SJed Brown dm->ltogmapb = dm->ltogmap; 2078ac4e037SJed Brown PetscFunctionReturn(0); 2088ac4e037SJed Brown } 2098ac4e037SJed Brown 2108ac4e037SJed Brown #undef __FUNCT__ 2118ac4e037SJed Brown #define __FUNCT__ "DMView_Redundant" 2128ac4e037SJed Brown static PetscErrorCode DMView_Redundant(DM dm,PetscViewer viewer) 2138ac4e037SJed Brown { 2148ac4e037SJed Brown PetscErrorCode ierr; 2158ac4e037SJed Brown DM_Redundant *red = (DM_Redundant*)dm->data; 2168ac4e037SJed Brown PetscBool iascii; 2178ac4e037SJed Brown 2188ac4e037SJed Brown PetscFunctionBegin; 219251f4c67SDmitry Karpeev ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 2208ac4e037SJed Brown if (iascii) { 2218ac4e037SJed Brown ierr = PetscViewerASCIIPrintf(viewer,"redundant: rank=%D N=%D\n",red->rank,red->N);CHKERRQ(ierr); 2228ac4e037SJed Brown } 2238ac4e037SJed Brown PetscFunctionReturn(0); 2248ac4e037SJed Brown } 2258ac4e037SJed Brown 22669527181SJed Brown #undef __FUNCT__ 227e727c939SJed Brown #define __FUNCT__ "DMCreateColoring_Redundant" 22819fd82e9SBarry Smith static PetscErrorCode DMCreateColoring_Redundant(DM dm,ISColoringType ctype,MatType mtype,ISColoring *coloring) 22969527181SJed Brown { 23069527181SJed Brown DM_Redundant *red = (DM_Redundant*)dm->data; 23169527181SJed Brown PetscErrorCode ierr; 23269527181SJed Brown PetscInt i,nloc; 23369527181SJed Brown ISColoringValue *colors; 23469527181SJed Brown 23569527181SJed Brown PetscFunctionBegin; 23669527181SJed Brown switch (ctype) { 23769527181SJed Brown case IS_COLORING_GLOBAL: 23869527181SJed Brown nloc = red->n; 23969527181SJed Brown break; 24069527181SJed Brown case IS_COLORING_GHOSTED: 24169527181SJed Brown nloc = red->N; 24269527181SJed Brown break; 243ce94432eSBarry Smith default: SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONG,"Unknown ISColoringType %d",(int)ctype); 24469527181SJed Brown } 24569527181SJed Brown ierr = PetscMalloc(nloc*sizeof(ISColoringValue),&colors);CHKERRQ(ierr); 24669527181SJed Brown for (i=0; i<nloc; i++) colors[i] = i; 247ce94432eSBarry Smith ierr = ISColoringCreate(PetscObjectComm((PetscObject)dm),red->N,nloc,colors,coloring);CHKERRQ(ierr); 24869527181SJed Brown ierr = ISColoringSetType(*coloring,ctype);CHKERRQ(ierr); 24969527181SJed Brown PetscFunctionReturn(0); 25069527181SJed Brown } 2518ac4e037SJed Brown 2528ac4e037SJed Brown #undef __FUNCT__ 2538ac4e037SJed Brown #define __FUNCT__ "DMRefine_Redundant" 2548ac4e037SJed Brown static PetscErrorCode DMRefine_Redundant(DM dmc,MPI_Comm comm,DM *dmf) 2558ac4e037SJed Brown { 2568ac4e037SJed Brown PetscErrorCode ierr; 2578ac4e037SJed Brown PetscMPIInt flag; 2588ac4e037SJed Brown DM_Redundant *redc = (DM_Redundant*)dmc->data; 2598ac4e037SJed Brown 2608ac4e037SJed Brown PetscFunctionBegin; 261ce94432eSBarry Smith if (comm == MPI_COMM_NULL) { 262ce94432eSBarry Smith ierr = PetscObjectGetComm((PetscObject)dmc,&comm);CHKERRQ(ierr); 263ce94432eSBarry Smith } 264ce94432eSBarry Smith ierr = MPI_Comm_compare(PetscObjectComm((PetscObject)dmc),comm,&flag);CHKERRQ(ierr); 265ce94432eSBarry Smith if (flag != MPI_CONGRUENT && flag != MPI_IDENT) SETERRQ(PetscObjectComm((PetscObject)dmc),PETSC_ERR_SUP,"cannot change communicators"); 2668ac4e037SJed Brown ierr = DMRedundantCreate(comm,redc->rank,redc->N,dmf);CHKERRQ(ierr); 2678ac4e037SJed Brown PetscFunctionReturn(0); 2688ac4e037SJed Brown } 2698ac4e037SJed Brown 2708ac4e037SJed Brown #undef __FUNCT__ 2718ac4e037SJed Brown #define __FUNCT__ "DMCoarsen_Redundant" 2728ac4e037SJed Brown static PetscErrorCode DMCoarsen_Redundant(DM dmf,MPI_Comm comm,DM *dmc) 2738ac4e037SJed Brown { 2748ac4e037SJed Brown PetscErrorCode ierr; 2758ac4e037SJed Brown PetscMPIInt flag; 2768ac4e037SJed Brown DM_Redundant *redf = (DM_Redundant*)dmf->data; 2778ac4e037SJed Brown 2788ac4e037SJed Brown PetscFunctionBegin; 279ce94432eSBarry Smith if (comm == MPI_COMM_NULL) { 280ce94432eSBarry Smith ierr = PetscObjectGetComm((PetscObject)dmf,&comm);CHKERRQ(ierr); 281ce94432eSBarry Smith } 282ce94432eSBarry Smith ierr = MPI_Comm_compare(PetscObjectComm((PetscObject)dmf),comm,&flag);CHKERRQ(ierr); 283ce94432eSBarry Smith if (flag != MPI_CONGRUENT && flag != MPI_IDENT) SETERRQ(PetscObjectComm((PetscObject)dmf),PETSC_ERR_SUP,"cannot change communicators"); 2848ac4e037SJed Brown ierr = DMRedundantCreate(comm,redf->rank,redf->N,dmc);CHKERRQ(ierr); 2858ac4e037SJed Brown PetscFunctionReturn(0); 2868ac4e037SJed Brown } 2878ac4e037SJed Brown 2888ac4e037SJed Brown #undef __FUNCT__ 289e727c939SJed Brown #define __FUNCT__ "DMCreateInterpolation_Redundant" 290e727c939SJed Brown static PetscErrorCode DMCreateInterpolation_Redundant(DM dmc,DM dmf,Mat *P,Vec *scale) 2918ac4e037SJed Brown { 2928ac4e037SJed Brown PetscErrorCode ierr; 2938ac4e037SJed Brown DM_Redundant *redc = (DM_Redundant*)dmc->data; 2948ac4e037SJed Brown DM_Redundant *redf = (DM_Redundant*)dmf->data; 2958ac4e037SJed Brown PetscMPIInt flag; 2968ac4e037SJed Brown PetscInt i,rstart,rend; 2978ac4e037SJed Brown 2988ac4e037SJed Brown PetscFunctionBegin; 299ce94432eSBarry Smith ierr = MPI_Comm_compare(PetscObjectComm((PetscObject)dmc),PetscObjectComm((PetscObject)dmf),&flag);CHKERRQ(ierr); 300ce94432eSBarry Smith if (flag != MPI_CONGRUENT && flag != MPI_IDENT) SETERRQ(PetscObjectComm((PetscObject)dmf),PETSC_ERR_SUP,"cannot change communicators"); 301ce94432eSBarry Smith if (redc->rank != redf->rank) SETERRQ(PetscObjectComm((PetscObject)dmf),PETSC_ERR_ARG_INCOMP,"Owning rank does not match"); 302ce94432eSBarry Smith if (redc->N != redf->N) SETERRQ(PetscObjectComm((PetscObject)dmf),PETSC_ERR_ARG_INCOMP,"Global size does not match"); 303ce94432eSBarry Smith ierr = MatCreate(PetscObjectComm((PetscObject)dmc),P);CHKERRQ(ierr); 3048ac4e037SJed Brown ierr = MatSetSizes(*P,redc->n,redc->n,redc->N,redc->N);CHKERRQ(ierr); 3058ac4e037SJed Brown ierr = MatSetType(*P,MATAIJ);CHKERRQ(ierr); 3068ac4e037SJed Brown ierr = MatSeqAIJSetPreallocation(*P,1,0);CHKERRQ(ierr); 3078ac4e037SJed Brown ierr = MatMPIAIJSetPreallocation(*P,1,0,0,0);CHKERRQ(ierr); 3088ac4e037SJed Brown ierr = MatGetOwnershipRange(*P,&rstart,&rend);CHKERRQ(ierr); 3098ac4e037SJed Brown for (i=rstart; i<rend; i++) {ierr = MatSetValue(*P,i,i,1.0,INSERT_VALUES);CHKERRQ(ierr);} 3108ac4e037SJed Brown ierr = MatAssemblyBegin(*P,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3118ac4e037SJed Brown ierr = MatAssemblyEnd(*P,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 312e727c939SJed Brown if (scale) {ierr = DMCreateInterpolationScale(dmc,dmf,*P,scale);CHKERRQ(ierr);} 3138ac4e037SJed Brown PetscFunctionReturn(0); 3148ac4e037SJed Brown } 3158ac4e037SJed Brown 3168ac4e037SJed Brown #undef __FUNCT__ 3178ac4e037SJed Brown #define __FUNCT__ "DMRedundantSetSize" 3188ac4e037SJed Brown /*@ 3198ac4e037SJed Brown DMRedundantSetSize - Sets the size of a densely coupled redundant object 3208ac4e037SJed Brown 3218ac4e037SJed Brown Collective on DM 3228ac4e037SJed Brown 3238ac4e037SJed Brown Input Parameter: 3248ac4e037SJed Brown + dm - redundant DM 3258ac4e037SJed Brown . rank - rank of process to own redundant degrees of freedom 3268ac4e037SJed Brown - N - total number of redundant degrees of freedom 3278ac4e037SJed Brown 3288ac4e037SJed Brown Level: advanced 3298ac4e037SJed Brown 3308ac4e037SJed Brown .seealso DMDestroy(), DMCreateGlobalVector(), DMRedundantCreate(), DMRedundantGetSize() 3318ac4e037SJed Brown @*/ 3328ac4e037SJed Brown PetscErrorCode DMRedundantSetSize(DM dm,PetscInt rank,PetscInt N) 3338ac4e037SJed Brown { 3348ac4e037SJed Brown PetscErrorCode ierr; 3358ac4e037SJed Brown 3368ac4e037SJed Brown PetscFunctionBegin; 3378ac4e037SJed Brown PetscValidHeaderSpecific(dm,DM_CLASSID,1); 3388ac4e037SJed Brown PetscValidType(dm,1); 3398ac4e037SJed Brown PetscValidLogicalCollectiveInt(dm,rank,2); 3408ac4e037SJed Brown PetscValidLogicalCollectiveInt(dm,N,3); 3418ac4e037SJed Brown ierr = PetscTryMethod(dm,"DMRedundantSetSize_C",(DM,PetscInt,PetscInt),(dm,rank,N));CHKERRQ(ierr); 3428ac4e037SJed Brown PetscFunctionReturn(0); 3438ac4e037SJed Brown } 3448ac4e037SJed Brown 3458ac4e037SJed Brown #undef __FUNCT__ 3468ac4e037SJed Brown #define __FUNCT__ "DMRedundantGetSize" 3478ac4e037SJed Brown /*@ 3488ac4e037SJed Brown DMRedundantGetSize - Gets the size of a densely coupled redundant object 3498ac4e037SJed Brown 3508ac4e037SJed Brown Not Collective 3518ac4e037SJed Brown 3528ac4e037SJed Brown Input Parameter: 3538ac4e037SJed Brown + dm - redundant DM 3548ac4e037SJed Brown 3558ac4e037SJed Brown Output Parameters: 3560298fd71SBarry Smith + rank - rank of process to own redundant degrees of freedom (or NULL) 3570298fd71SBarry Smith - N - total number of redundant degrees of freedom (or NULL) 3588ac4e037SJed Brown 3598ac4e037SJed Brown Level: advanced 3608ac4e037SJed Brown 3618ac4e037SJed Brown .seealso DMDestroy(), DMCreateGlobalVector(), DMRedundantCreate(), DMRedundantSetSize() 3628ac4e037SJed Brown @*/ 3638ac4e037SJed Brown PetscErrorCode DMRedundantGetSize(DM dm,PetscInt *rank,PetscInt *N) 3648ac4e037SJed Brown { 3658ac4e037SJed Brown PetscErrorCode ierr; 3668ac4e037SJed Brown 3678ac4e037SJed Brown PetscFunctionBegin; 3688ac4e037SJed Brown PetscValidHeaderSpecific(dm,DM_CLASSID,1); 3698ac4e037SJed Brown PetscValidType(dm,1); 3708ac4e037SJed Brown ierr = PetscUseMethod(dm,"DMRedundantGetSize_C",(DM,PetscInt*,PetscInt*),(dm,rank,N));CHKERRQ(ierr); 3718ac4e037SJed Brown PetscFunctionReturn(0); 3728ac4e037SJed Brown } 3738ac4e037SJed Brown 3748ac4e037SJed Brown #undef __FUNCT__ 3758ac4e037SJed Brown #define __FUNCT__ "DMRedundantSetSize_Redundant" 3761e6b0712SBarry Smith static PetscErrorCode DMRedundantSetSize_Redundant(DM dm,PetscInt rank,PetscInt N) 3778ac4e037SJed Brown { 3788ac4e037SJed Brown DM_Redundant *red = (DM_Redundant*)dm->data; 3798ac4e037SJed Brown PetscErrorCode ierr; 3808ac4e037SJed Brown PetscMPIInt myrank; 3818ac4e037SJed Brown 3828ac4e037SJed Brown PetscFunctionBegin; 383ce94432eSBarry Smith ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&myrank);CHKERRQ(ierr); 3848ac4e037SJed Brown red->rank = rank; 3858ac4e037SJed Brown red->N = N; 3868ac4e037SJed Brown red->n = (myrank == rank) ? N : 0; 3878ac4e037SJed Brown PetscFunctionReturn(0); 3888ac4e037SJed Brown } 3898ac4e037SJed Brown 3908ac4e037SJed Brown #undef __FUNCT__ 3918ac4e037SJed Brown #define __FUNCT__ "DMRedundantGetSize_Redundant" 3921e6b0712SBarry Smith static PetscErrorCode DMRedundantGetSize_Redundant(DM dm,PetscInt *rank,PetscInt *N) 3938ac4e037SJed Brown { 3948ac4e037SJed Brown DM_Redundant *red = (DM_Redundant*)dm->data; 3958ac4e037SJed Brown 3968ac4e037SJed Brown PetscFunctionBegin; 3978ac4e037SJed Brown if (rank) *rank = red->rank; 3988ac4e037SJed Brown if (N) *N = red->N; 3998ac4e037SJed Brown PetscFunctionReturn(0); 4008ac4e037SJed Brown } 4018ac4e037SJed Brown 4023efe6655SBarry Smith /*MC 4033efe6655SBarry Smith DMREDUNDANT = "redundant" - A DM object that is used to manage data for a small set of dense globally coupled variables. 4043efe6655SBarry Smith In the global representation of the vector the variables are all stored on a single MPI process (all the other MPI processes have 4053efe6655SBarry Smith no variables) in the local representation all the variables are stored on ALL the MPI processes (because they are all needed for each 4063efe6655SBarry Smith processes local computations). 4073efe6655SBarry Smith 4083efe6655SBarry Smith This DM is generally used inside a DMCOMPOSITE object. For example, it may be used to store continuation parameters for a bifurcation problem. 4093efe6655SBarry Smith 4103efe6655SBarry Smith 4113efe6655SBarry Smith Level: intermediate 4123efe6655SBarry Smith 4133efe6655SBarry Smith .seealso: DMType, DMCOMPOSITE, DMCreateRedundant(), DMCreate(), DMRedundantSetSize(), DMRedundantGetSize() 4143efe6655SBarry Smith M*/ 4153efe6655SBarry Smith 4168ac4e037SJed Brown #undef __FUNCT__ 4178ac4e037SJed Brown #define __FUNCT__ "DMCreate_Redundant" 418*8cc058d9SJed Brown PETSC_EXTERN PetscErrorCode DMCreate_Redundant(DM dm) 4198ac4e037SJed Brown { 4208ac4e037SJed Brown PetscErrorCode ierr; 4218ac4e037SJed Brown DM_Redundant *red; 4228ac4e037SJed Brown 4238ac4e037SJed Brown PetscFunctionBegin; 4248ac4e037SJed Brown ierr = PetscNewLog(dm,DM_Redundant,&red);CHKERRQ(ierr); 4258ac4e037SJed Brown dm->data = red; 4268ac4e037SJed Brown 4278ac4e037SJed Brown ierr = PetscObjectChangeTypeName((PetscObject)dm,DMREDUNDANT);CHKERRQ(ierr); 4288865f1eaSKarl Rupp 4298ac4e037SJed Brown dm->ops->setup = DMSetUp_Redundant; 4308ac4e037SJed Brown dm->ops->view = DMView_Redundant; 4318ac4e037SJed Brown dm->ops->createglobalvector = DMCreateGlobalVector_Redundant; 4328ac4e037SJed Brown dm->ops->createlocalvector = DMCreateLocalVector_Redundant; 43325296bd5SBarry Smith dm->ops->creatematrix = DMCreateMatrix_Redundant; 4348ac4e037SJed Brown dm->ops->destroy = DMDestroy_Redundant; 4358ac4e037SJed Brown dm->ops->globaltolocalbegin = DMGlobalToLocalBegin_Redundant; 4368ac4e037SJed Brown dm->ops->globaltolocalend = DMGlobalToLocalEnd_Redundant; 4378ac4e037SJed Brown dm->ops->localtoglobalbegin = DMLocalToGlobalBegin_Redundant; 4388ac4e037SJed Brown dm->ops->localtoglobalend = DMLocalToGlobalEnd_Redundant; 4398ac4e037SJed Brown dm->ops->refine = DMRefine_Redundant; 4408ac4e037SJed Brown dm->ops->coarsen = DMCoarsen_Redundant; 44125296bd5SBarry Smith dm->ops->createinterpolation= DMCreateInterpolation_Redundant; 442e727c939SJed Brown dm->ops->getcoloring = DMCreateColoring_Redundant; 4438865f1eaSKarl Rupp 44419fd82e9SBarry Smith ierr = PetscStrallocpy(VECSTANDARD,(char**)&dm->vectype);CHKERRQ(ierr); 44500de8ff0SBarry Smith ierr = PetscObjectComposeFunction((PetscObject)dm,"DMRedundantSetSize_C","DMRedundantSetSize_Redundant",DMRedundantSetSize_Redundant);CHKERRQ(ierr); 44600de8ff0SBarry Smith ierr = PetscObjectComposeFunction((PetscObject)dm,"DMRedundantGetSize_C","DMRedundantGetSize_Redundant",DMRedundantGetSize_Redundant);CHKERRQ(ierr); 4478ac4e037SJed Brown PetscFunctionReturn(0); 4488ac4e037SJed Brown } 4498ac4e037SJed Brown 4508ac4e037SJed Brown #undef __FUNCT__ 4518ac4e037SJed Brown #define __FUNCT__ "DMRedundantCreate" 4528ac4e037SJed Brown /*@C 4538ac4e037SJed Brown DMRedundantCreate - Creates a DM object, used to manage data for dense globally coupled variables 4548ac4e037SJed Brown 4558ac4e037SJed Brown Collective on MPI_Comm 4568ac4e037SJed Brown 4578ac4e037SJed Brown Input Parameter: 4588ac4e037SJed Brown + comm - the processors that will share the global vector 4598ac4e037SJed Brown . rank - rank to own the redundant values 4608ac4e037SJed Brown - N - total number of degrees of freedom 4618ac4e037SJed Brown 4628ac4e037SJed Brown Output Parameters: 4638ac4e037SJed Brown . red - the redundant DM 4648ac4e037SJed Brown 4658ac4e037SJed Brown Level: advanced 4668ac4e037SJed Brown 4673efe6655SBarry Smith .seealso DMDestroy(), DMCreateGlobalVector(), DMCreateMatrix(), DMCompositeAddDM(), DMREDUNDANT, DMSetType(), DMRedundantSetSize(), DMRedundantGetSize() 4688ac4e037SJed Brown 4698ac4e037SJed Brown @*/ 4708ac4e037SJed Brown PetscErrorCode DMRedundantCreate(MPI_Comm comm,PetscInt rank,PetscInt N,DM *dm) 4718ac4e037SJed Brown { 4728ac4e037SJed Brown PetscErrorCode ierr; 4738ac4e037SJed Brown 4748ac4e037SJed Brown PetscFunctionBegin; 4758ac4e037SJed Brown PetscValidPointer(dm,2); 4768ac4e037SJed Brown ierr = DMCreate(comm,dm);CHKERRQ(ierr); 4778ac4e037SJed Brown ierr = DMSetType(*dm,DMREDUNDANT);CHKERRQ(ierr); 4788ac4e037SJed Brown ierr = DMRedundantSetSize(*dm,rank,N);CHKERRQ(ierr); 4798ac4e037SJed Brown ierr = DMSetUp(*dm);CHKERRQ(ierr); 4808ac4e037SJed Brown PetscFunctionReturn(0); 4818ac4e037SJed Brown } 482