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