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