1 2 #include <../src/dm/impls/composite/packimpl.h> /*I "petscdmcomposite.h" I*/ 3 4 #undef __FUNCT__ 5 #define __FUNCT__ "DMCreateMatrix_Composite_Nest" 6 static PetscErrorCode DMCreateMatrix_Composite_Nest(DM dm,MatType mtype,Mat *J) 7 { 8 const DM_Composite *com = (DM_Composite*)dm->data; 9 const struct DMCompositeLink *rlink,*clink; 10 PetscErrorCode ierr; 11 IS *isg; 12 Mat *submats; 13 PetscInt i,j,n; 14 15 PetscFunctionBegin; 16 n = com->nDM; /* Total number of entries */ 17 18 /* Explicit index sets are not required for MatCreateNest, but getting them here allows MatNest to do consistency 19 * checking and allows ISEqual to compare by identity instead of by contents. */ 20 ierr = DMCompositeGetGlobalISs(dm,&isg);CHKERRQ(ierr); 21 22 /* Get submatrices */ 23 ierr = PetscMalloc(n*n*sizeof(Mat),&submats);CHKERRQ(ierr); 24 for (i=0,rlink=com->next; rlink; i++,rlink=rlink->next) { 25 for (j=0,clink=com->next; clink; j++,clink=clink->next) { 26 Mat sub = PETSC_NULL; 27 if (i == j) { 28 ierr = DMCreateMatrix(rlink->dm,PETSC_NULL,&sub);CHKERRQ(ierr); 29 } else if (com->FormCoupleLocations) SETERRQ(((PetscObject)dm)->comm,PETSC_ERR_SUP,"Cannot manage off-diagonal parts yet"); 30 submats[i*n+j] = sub; 31 } 32 } 33 34 ierr = MatCreateNest(((PetscObject)dm)->comm,n,isg,n,isg,submats,J);CHKERRQ(ierr); 35 36 /* Disown references */ 37 for (i=0; i<n; i++) {ierr = ISDestroy(&isg[i]);CHKERRQ(ierr);} 38 ierr = PetscFree(isg);CHKERRQ(ierr); 39 40 for (i=0; i<n*n; i++) { 41 if (submats[i]) {ierr = MatDestroy(&submats[i]);CHKERRQ(ierr);} 42 } 43 ierr = PetscFree(submats);CHKERRQ(ierr); 44 PetscFunctionReturn(0); 45 } 46 47 #undef __FUNCT__ 48 #define __FUNCT__ "DMCreateMatrix_Composite_AIJ" 49 static PetscErrorCode DMCreateMatrix_Composite_AIJ(DM dm,MatType mtype,Mat *J) 50 { 51 PetscErrorCode ierr; 52 DM_Composite *com = (DM_Composite*)dm->data; 53 struct DMCompositeLink *next; 54 PetscInt m,*dnz,*onz,i,j,mA; 55 Mat Atmp; 56 PetscMPIInt rank; 57 PetscBool dense = PETSC_FALSE; 58 59 PetscFunctionBegin; 60 /* use global vector to determine layout needed for matrix */ 61 m = com->n; 62 63 ierr = MatCreate(((PetscObject)dm)->comm,J);CHKERRQ(ierr); 64 ierr = MatSetSizes(*J,m,m,PETSC_DETERMINE,PETSC_DETERMINE);CHKERRQ(ierr); 65 ierr = MatSetType(*J,mtype);CHKERRQ(ierr); 66 67 /* 68 Extremely inefficient but will compute entire Jacobian for testing 69 */ 70 ierr = PetscOptionsGetBool(((PetscObject)dm)->prefix,"-dmcomposite_dense_jacobian",&dense,PETSC_NULL);CHKERRQ(ierr); 71 if (dense) { 72 PetscInt rstart,rend,*indices; 73 PetscScalar *values; 74 75 mA = com->N; 76 ierr = MatMPIAIJSetPreallocation(*J,mA,PETSC_NULL,mA-m,PETSC_NULL);CHKERRQ(ierr); 77 ierr = MatSeqAIJSetPreallocation(*J,mA,PETSC_NULL);CHKERRQ(ierr); 78 79 ierr = MatGetOwnershipRange(*J,&rstart,&rend);CHKERRQ(ierr); 80 ierr = PetscMalloc2(mA,PetscScalar,&values,mA,PetscInt,&indices);CHKERRQ(ierr); 81 ierr = PetscMemzero(values,mA*sizeof(PetscScalar));CHKERRQ(ierr); 82 for (i=0; i<mA; i++) indices[i] = i; 83 for (i=rstart; i<rend; i++) { 84 ierr = MatSetValues(*J,1,&i,mA,indices,values,INSERT_VALUES);CHKERRQ(ierr); 85 } 86 ierr = PetscFree2(values,indices);CHKERRQ(ierr); 87 ierr = MatAssemblyBegin(*J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 88 ierr = MatAssemblyEnd(*J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 89 PetscFunctionReturn(0); 90 } 91 92 ierr = MPI_Comm_rank(((PetscObject)dm)->comm,&rank);CHKERRQ(ierr); 93 ierr = MatPreallocateInitialize(((PetscObject)dm)->comm,m,m,dnz,onz);CHKERRQ(ierr); 94 /* loop over packed objects, handling one at at time */ 95 next = com->next; 96 while (next) { 97 PetscInt nc,rstart,*ccols,maxnc; 98 const PetscInt *cols,*rstarts; 99 PetscMPIInt proc; 100 101 ierr = DMCreateMatrix(next->dm,mtype,&Atmp);CHKERRQ(ierr); 102 ierr = MatGetOwnershipRange(Atmp,&rstart,PETSC_NULL);CHKERRQ(ierr); 103 ierr = MatGetOwnershipRanges(Atmp,&rstarts);CHKERRQ(ierr); 104 ierr = MatGetLocalSize(Atmp,&mA,PETSC_NULL);CHKERRQ(ierr); 105 106 maxnc = 0; 107 for (i=0; i<mA; i++) { 108 ierr = MatGetRow(Atmp,rstart+i,&nc,PETSC_NULL,PETSC_NULL);CHKERRQ(ierr); 109 ierr = MatRestoreRow(Atmp,rstart+i,&nc,PETSC_NULL,PETSC_NULL);CHKERRQ(ierr); 110 maxnc = PetscMax(nc,maxnc); 111 } 112 ierr = PetscMalloc(maxnc*sizeof(PetscInt),&ccols);CHKERRQ(ierr); 113 for (i=0; i<mA; i++) { 114 ierr = MatGetRow(Atmp,rstart+i,&nc,&cols,PETSC_NULL);CHKERRQ(ierr); 115 /* remap the columns taking into how much they are shifted on each process */ 116 for (j=0; j<nc; j++) { 117 proc = 0; 118 while (cols[j] >= rstarts[proc+1]) proc++; 119 ccols[j] = cols[j] + next->grstarts[proc] - rstarts[proc]; 120 } 121 ierr = MatPreallocateSet(com->rstart+next->rstart+i,nc,ccols,dnz,onz);CHKERRQ(ierr); 122 ierr = MatRestoreRow(Atmp,rstart+i,&nc,&cols,PETSC_NULL);CHKERRQ(ierr); 123 } 124 ierr = PetscFree(ccols);CHKERRQ(ierr); 125 ierr = MatDestroy(&Atmp);CHKERRQ(ierr); 126 next = next->next; 127 } 128 if (com->FormCoupleLocations) { 129 ierr = (*com->FormCoupleLocations)(dm,PETSC_NULL,dnz,onz,__rstart,__nrows,__start,__end);CHKERRQ(ierr); 130 } 131 ierr = MatMPIAIJSetPreallocation(*J,0,dnz,0,onz);CHKERRQ(ierr); 132 ierr = MatSeqAIJSetPreallocation(*J,0,dnz);CHKERRQ(ierr); 133 ierr = MatPreallocateFinalize(dnz,onz);CHKERRQ(ierr); 134 135 if (dm->prealloc_only) PetscFunctionReturn(0); 136 137 next = com->next; 138 while (next) { 139 PetscInt nc,rstart,row,maxnc,*ccols; 140 const PetscInt *cols,*rstarts; 141 const PetscScalar *values; 142 PetscMPIInt proc; 143 144 ierr = DMCreateMatrix(next->dm,mtype,&Atmp);CHKERRQ(ierr); 145 ierr = MatGetOwnershipRange(Atmp,&rstart,PETSC_NULL);CHKERRQ(ierr); 146 ierr = MatGetOwnershipRanges(Atmp,&rstarts);CHKERRQ(ierr); 147 ierr = MatGetLocalSize(Atmp,&mA,PETSC_NULL);CHKERRQ(ierr); 148 maxnc = 0; 149 for (i=0; i<mA; i++) { 150 ierr = MatGetRow(Atmp,rstart+i,&nc,PETSC_NULL,PETSC_NULL);CHKERRQ(ierr); 151 ierr = MatRestoreRow(Atmp,rstart+i,&nc,PETSC_NULL,PETSC_NULL);CHKERRQ(ierr); 152 maxnc = PetscMax(nc,maxnc); 153 } 154 ierr = PetscMalloc(maxnc*sizeof(PetscInt),&ccols);CHKERRQ(ierr); 155 for (i=0; i<mA; i++) { 156 ierr = MatGetRow(Atmp,rstart+i,&nc,(const PetscInt**)&cols,&values);CHKERRQ(ierr); 157 for (j=0; j<nc; j++) { 158 proc = 0; 159 while (cols[j] >= rstarts[proc+1]) proc++; 160 ccols[j] = cols[j] + next->grstarts[proc] - rstarts[proc]; 161 } 162 row = com->rstart+next->rstart+i; 163 ierr = MatSetValues(*J,1,&row,nc,ccols,values,INSERT_VALUES);CHKERRQ(ierr); 164 ierr = MatRestoreRow(Atmp,rstart+i,&nc,(const PetscInt**)&cols,&values);CHKERRQ(ierr); 165 } 166 ierr = PetscFree(ccols);CHKERRQ(ierr); 167 ierr = MatDestroy(&Atmp);CHKERRQ(ierr); 168 next = next->next; 169 } 170 if (com->FormCoupleLocations) { 171 PetscInt __rstart; 172 ierr = MatGetOwnershipRange(*J,&__rstart,PETSC_NULL);CHKERRQ(ierr); 173 ierr = (*com->FormCoupleLocations)(dm,*J,PETSC_NULL,PETSC_NULL,__rstart,0,0,0);CHKERRQ(ierr); 174 } 175 ierr = MatAssemblyBegin(*J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 176 ierr = MatAssemblyEnd(*J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 177 PetscFunctionReturn(0); 178 } 179 180 #undef __FUNCT__ 181 #define __FUNCT__ "DMCreateMatrix_Composite" 182 PetscErrorCode DMCreateMatrix_Composite(DM dm,MatType mtype,Mat *J) 183 { 184 PetscErrorCode ierr; 185 PetscBool usenest; 186 ISLocalToGlobalMapping ltogmap,ltogmapb; 187 188 PetscFunctionBegin; 189 ierr = DMSetUp(dm);CHKERRQ(ierr); 190 ierr = PetscStrcmp(mtype,MATNEST,&usenest);CHKERRQ(ierr); 191 if (usenest) { 192 ierr = DMCreateMatrix_Composite_Nest(dm,mtype,J);CHKERRQ(ierr); 193 } else { 194 ierr = DMCreateMatrix_Composite_AIJ(dm,mtype ? mtype : MATAIJ,J);CHKERRQ(ierr); 195 } 196 197 ierr = DMGetLocalToGlobalMapping(dm,<ogmap);CHKERRQ(ierr); 198 ierr = DMGetLocalToGlobalMappingBlock(dm,<ogmapb);CHKERRQ(ierr); 199 ierr = MatSetLocalToGlobalMapping(*J,ltogmap,ltogmap);CHKERRQ(ierr); 200 ierr = MatSetLocalToGlobalMappingBlock(*J,ltogmapb,ltogmapb);CHKERRQ(ierr); 201 PetscFunctionReturn(0); 202 } 203