1 2 #include "packimpl.h" /*I "petscdmcomposite.h" I*/ 3 4 #undef __FUNCT__ 5 #define __FUNCT__ "DMGetMatrix_Composite_Nest" 6 static PetscErrorCode DMGetMatrix_Composite_Nest(DM dm,const 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 + com->nredundant; /* 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 switch (rlink->type) { 29 case DMCOMPOSITE_ARRAY: 30 ierr = MatCreateMPIDense(((PetscObject)dm)->comm,rlink->n,clink->n,PETSC_DETERMINE,PETSC_DETERMINE,PETSC_NULL,&sub);CHKERRQ(ierr); 31 break; 32 case DMCOMPOSITE_DM: 33 ierr = DMGetMatrix(rlink->dm,PETSC_NULL,&sub);CHKERRQ(ierr); 34 break; 35 default: SETERRQ(((PetscObject)dm)->comm,PETSC_ERR_SUP,"unknown type"); 36 } 37 } else if (com->FormCoupleLocations) { 38 SETERRQ(((PetscObject)dm)->comm,PETSC_ERR_SUP,"Cannot manage off-diagonal parts yet"); 39 } 40 submats[i*n+j] = sub; 41 } 42 } 43 44 ierr = MatCreateNest(((PetscObject)dm)->comm,n,isg,n,isg,submats,J);CHKERRQ(ierr); 45 46 /* Disown references */ 47 for (i=0; i<n; i++) {ierr = ISDestroy(&isg[i]);CHKERRQ(ierr);} 48 ierr = PetscFree(isg);CHKERRQ(ierr); 49 50 for (i=0; i<n*n; i++) { 51 if (submats[i]) {ierr = MatDestroy(&submats[i]);CHKERRQ(ierr);} 52 } 53 ierr = PetscFree(submats);CHKERRQ(ierr); 54 PetscFunctionReturn(0); 55 } 56 57 #undef __FUNCT__ 58 #define __FUNCT__ "DMGetMatrix_Composite_AIJ" 59 static PetscErrorCode DMGetMatrix_Composite_AIJ(DM dm,const MatType mtype,Mat *J) 60 { 61 PetscErrorCode ierr; 62 DM_Composite *com = (DM_Composite*)dm->data; 63 struct DMCompositeLink *next = com->next; 64 PetscInt m,*dnz,*onz,i,j,mA; 65 Mat Atmp; 66 PetscMPIInt rank; 67 PetscBool dense = PETSC_FALSE; 68 69 PetscFunctionBegin; 70 /* use global vector to determine layout needed for matrix */ 71 m = com->n; 72 73 ierr = MatCreate(((PetscObject)dm)->comm,J);CHKERRQ(ierr); 74 ierr = MatSetSizes(*J,m,m,PETSC_DETERMINE,PETSC_DETERMINE);CHKERRQ(ierr); 75 ierr = MatSetType(*J,mtype);CHKERRQ(ierr); 76 77 /* 78 Extremely inefficient but will compute entire Jacobian for testing 79 */ 80 ierr = PetscOptionsGetBool(((PetscObject)dm)->prefix,"-dmcomposite_dense_jacobian",&dense,PETSC_NULL);CHKERRQ(ierr); 81 if (dense) { 82 PetscInt rstart,rend,*indices; 83 PetscScalar *values; 84 85 mA = com->N; 86 ierr = MatMPIAIJSetPreallocation(*J,mA,PETSC_NULL,mA-m,PETSC_NULL);CHKERRQ(ierr); 87 ierr = MatSeqAIJSetPreallocation(*J,mA,PETSC_NULL);CHKERRQ(ierr); 88 89 ierr = MatGetOwnershipRange(*J,&rstart,&rend);CHKERRQ(ierr); 90 ierr = PetscMalloc2(mA,PetscScalar,&values,mA,PetscInt,&indices);CHKERRQ(ierr); 91 ierr = PetscMemzero(values,mA*sizeof(PetscScalar));CHKERRQ(ierr); 92 for (i=0; i<mA; i++) indices[i] = i; 93 for (i=rstart; i<rend; i++) { 94 ierr = MatSetValues(*J,1,&i,mA,indices,values,INSERT_VALUES);CHKERRQ(ierr); 95 } 96 ierr = PetscFree2(values,indices);CHKERRQ(ierr); 97 ierr = MatAssemblyBegin(*J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 98 ierr = MatAssemblyEnd(*J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 99 PetscFunctionReturn(0); 100 } 101 102 ierr = MPI_Comm_rank(((PetscObject)dm)->comm,&rank);CHKERRQ(ierr); 103 ierr = MatPreallocateInitialize(((PetscObject)dm)->comm,m,m,dnz,onz);CHKERRQ(ierr); 104 /* loop over packed objects, handling one at at time */ 105 next = com->next; 106 while (next) { 107 if (next->type == DMCOMPOSITE_ARRAY) { 108 if (rank == next->rank) { /* zero the "little" block */ 109 for (j=com->rstart+next->rstart; j<com->rstart+next->rstart+next->n; j++) { 110 for (i=com->rstart+next->rstart; i<com->rstart+next->rstart+next->n; i++) { 111 ierr = MatPreallocateSet(j,1,&i,dnz,onz);CHKERRQ(ierr); 112 } 113 } 114 } 115 } else if (next->type == DMCOMPOSITE_DM) { 116 PetscInt nc,rstart,*ccols,maxnc; 117 const PetscInt *cols,*rstarts; 118 PetscMPIInt proc; 119 120 ierr = DMGetMatrix(next->dm,mtype,&Atmp);CHKERRQ(ierr); 121 ierr = MatGetOwnershipRange(Atmp,&rstart,PETSC_NULL);CHKERRQ(ierr); 122 ierr = MatGetOwnershipRanges(Atmp,&rstarts);CHKERRQ(ierr); 123 ierr = MatGetLocalSize(Atmp,&mA,PETSC_NULL);CHKERRQ(ierr); 124 125 maxnc = 0; 126 for (i=0; i<mA; i++) { 127 ierr = MatGetRow(Atmp,rstart+i,&nc,PETSC_NULL,PETSC_NULL);CHKERRQ(ierr); 128 ierr = MatRestoreRow(Atmp,rstart+i,&nc,PETSC_NULL,PETSC_NULL);CHKERRQ(ierr); 129 maxnc = PetscMax(nc,maxnc); 130 } 131 ierr = PetscMalloc(maxnc*sizeof(PetscInt),&ccols);CHKERRQ(ierr); 132 for (i=0; i<mA; i++) { 133 ierr = MatGetRow(Atmp,rstart+i,&nc,&cols,PETSC_NULL);CHKERRQ(ierr); 134 /* remap the columns taking into how much they are shifted on each process */ 135 for (j=0; j<nc; j++) { 136 proc = 0; 137 while (cols[j] >= rstarts[proc+1]) proc++; 138 ccols[j] = cols[j] + next->grstarts[proc] - rstarts[proc]; 139 } 140 ierr = MatPreallocateSet(com->rstart+next->rstart+i,nc,ccols,dnz,onz);CHKERRQ(ierr); 141 ierr = MatRestoreRow(Atmp,rstart+i,&nc,&cols,PETSC_NULL);CHKERRQ(ierr); 142 } 143 ierr = PetscFree(ccols);CHKERRQ(ierr); 144 ierr = MatDestroy(&Atmp);CHKERRQ(ierr); 145 } else { 146 SETERRQ(((PetscObject)dm)->comm,PETSC_ERR_SUP,"Cannot handle that object type yet"); 147 } 148 next = next->next; 149 } 150 if (com->FormCoupleLocations) { 151 ierr = (*com->FormCoupleLocations)(dm,PETSC_NULL,dnz,onz,__rstart,__nrows,__start,__end);CHKERRQ(ierr); 152 } 153 ierr = MatMPIAIJSetPreallocation(*J,0,dnz,0,onz);CHKERRQ(ierr); 154 ierr = MatSeqAIJSetPreallocation(*J,0,dnz);CHKERRQ(ierr); 155 ierr = MatPreallocateFinalize(dnz,onz);CHKERRQ(ierr); 156 157 if (dm->prealloc_only) PetscFunctionReturn(0); 158 159 next = com->next; 160 while (next) { 161 if (next->type == DMCOMPOSITE_ARRAY) { 162 if (rank == next->rank) { 163 for (j=com->rstart+next->rstart; j<com->rstart+next->rstart+next->n; j++) { 164 for (i=com->rstart+next->rstart; i<com->rstart+next->rstart+next->n; i++) { 165 ierr = MatSetValue(*J,j,i,0.0,INSERT_VALUES);CHKERRQ(ierr); 166 } 167 } 168 } 169 } else if (next->type == DMCOMPOSITE_DM) { 170 PetscInt nc,rstart,row,maxnc,*ccols; 171 const PetscInt *cols,*rstarts; 172 const PetscScalar *values; 173 PetscMPIInt proc; 174 175 ierr = DMGetMatrix(next->dm,mtype,&Atmp);CHKERRQ(ierr); 176 ierr = MatGetOwnershipRange(Atmp,&rstart,PETSC_NULL);CHKERRQ(ierr); 177 ierr = MatGetOwnershipRanges(Atmp,&rstarts);CHKERRQ(ierr); 178 ierr = MatGetLocalSize(Atmp,&mA,PETSC_NULL);CHKERRQ(ierr); 179 maxnc = 0; 180 for (i=0; i<mA; i++) { 181 ierr = MatGetRow(Atmp,rstart+i,&nc,PETSC_NULL,PETSC_NULL);CHKERRQ(ierr); 182 ierr = MatRestoreRow(Atmp,rstart+i,&nc,PETSC_NULL,PETSC_NULL);CHKERRQ(ierr); 183 maxnc = PetscMax(nc,maxnc); 184 } 185 ierr = PetscMalloc(maxnc*sizeof(PetscInt),&ccols);CHKERRQ(ierr); 186 for (i=0; i<mA; i++) { 187 ierr = MatGetRow(Atmp,rstart+i,&nc,(const PetscInt **)&cols,&values);CHKERRQ(ierr); 188 for (j=0; j<nc; j++) { 189 proc = 0; 190 while (cols[j] >= rstarts[proc+1]) proc++; 191 ccols[j] = cols[j] + next->grstarts[proc] - rstarts[proc]; 192 } 193 row = com->rstart+next->rstart+i; 194 ierr = MatSetValues(*J,1,&row,nc,ccols,values,INSERT_VALUES);CHKERRQ(ierr); 195 ierr = MatRestoreRow(Atmp,rstart+i,&nc,(const PetscInt **)&cols,&values);CHKERRQ(ierr); 196 } 197 ierr = PetscFree(ccols);CHKERRQ(ierr); 198 ierr = MatDestroy(&Atmp);CHKERRQ(ierr); 199 } else { 200 SETERRQ(((PetscObject)dm)->comm,PETSC_ERR_SUP,"Cannot handle that object type yet"); 201 } 202 next = next->next; 203 } 204 if (com->FormCoupleLocations) { 205 PetscInt __rstart; 206 ierr = MatGetOwnershipRange(*J,&__rstart,PETSC_NULL);CHKERRQ(ierr); 207 ierr = (*com->FormCoupleLocations)(dm,*J,PETSC_NULL,PETSC_NULL,__rstart,0,0,0);CHKERRQ(ierr); 208 } 209 ierr = MatAssemblyBegin(*J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 210 ierr = MatAssemblyEnd(*J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 211 PetscFunctionReturn(0); 212 } 213 214 #undef __FUNCT__ 215 #define __FUNCT__ "DMGetMatrix_Composite" 216 PetscErrorCode DMGetMatrix_Composite(DM dm,const MatType mtype,Mat *J) 217 { 218 PetscErrorCode ierr; 219 PetscBool usenest; 220 ISLocalToGlobalMapping ltogmap,ltogmapb; 221 222 PetscFunctionBegin; 223 ierr = PetscStrcmp(mtype,MATNEST,&usenest);CHKERRQ(ierr); 224 if (usenest) { 225 ierr = DMGetMatrix_Composite_Nest(dm,mtype,J);CHKERRQ(ierr); 226 } else { 227 ierr = DMGetMatrix_Composite_AIJ(dm,mtype?mtype:MATAIJ,J);CHKERRQ(ierr); 228 } 229 230 ierr = DMGetLocalToGlobalMapping(dm,<ogmap);CHKERRQ(ierr); 231 ierr = DMGetLocalToGlobalMappingBlock(dm,<ogmapb);CHKERRQ(ierr); 232 ierr = MatSetLocalToGlobalMapping(*J,ltogmap,ltogmap);CHKERRQ(ierr); 233 ierr = MatSetLocalToGlobalMappingBlock(*J,ltogmapb,ltogmapb);CHKERRQ(ierr); 234 PetscFunctionReturn(0); 235 } 236