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