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