1 2 #include <petsc-private/matimpl.h> /*I "petscmat.h" I*/ 3 4 typedef struct { 5 Mat Top; 6 PetscErrorCode (*SetValues)(Mat,PetscInt,const PetscInt[],PetscInt,const PetscInt[],const PetscScalar[],InsertMode); 7 PetscErrorCode (*SetValuesBlocked)(Mat,PetscInt,const PetscInt[],PetscInt,const PetscInt[],const PetscScalar[],InsertMode); 8 } Mat_LocalRef; 9 10 /* These need to be macros because they use sizeof */ 11 #define IndexSpaceGet(buf,nrow,ncol,irowm,icolm) do { \ 12 if (nrow + ncol > (PetscInt)(sizeof(buf)/sizeof(buf[0]))) { \ 13 ierr = PetscMalloc2(nrow,&irowm,ncol,&icolm);CHKERRQ(ierr); \ 14 } else { \ 15 irowm = &buf[0]; \ 16 icolm = &buf[nrow]; \ 17 } \ 18 } while (0) 19 20 #define IndexSpaceRestore(buf,nrow,ncol,irowm,icolm) do { \ 21 if (nrow + ncol > (PetscInt)(sizeof(buf)/sizeof(buf[0]))) { \ 22 ierr = PetscFree2(irowm,icolm);CHKERRQ(ierr); \ 23 } \ 24 } while (0) 25 26 static void BlockIndicesExpand(PetscInt n,const PetscInt idx[],PetscInt bs,PetscInt idxm[]) 27 { 28 PetscInt i,j; 29 for (i=0; i<n; i++) { 30 for (j=0; j<bs; j++) { 31 idxm[i*bs+j] = idx[i]*bs + j; 32 } 33 } 34 } 35 36 #undef __FUNCT__ 37 #define __FUNCT__ "MatSetValuesBlockedLocal_LocalRef_Block" 38 static PetscErrorCode MatSetValuesBlockedLocal_LocalRef_Block(Mat A,PetscInt nrow,const PetscInt irow[],PetscInt ncol,const PetscInt icol[],const PetscScalar y[],InsertMode addv) 39 { 40 Mat_LocalRef *lr = (Mat_LocalRef*)A->data; 41 PetscErrorCode ierr; 42 PetscInt buf[4096],*irowm,*icolm; 43 44 PetscFunctionBegin; 45 if (!nrow || !ncol) PetscFunctionReturn(0); 46 IndexSpaceGet(buf,nrow,ncol,irowm,icolm); 47 ierr = ISLocalToGlobalMappingApplyBlock(A->rmap->mapping,nrow,irow,irowm);CHKERRQ(ierr); 48 ierr = ISLocalToGlobalMappingApplyBlock(A->cmap->mapping,ncol,icol,icolm);CHKERRQ(ierr); 49 ierr = (*lr->SetValuesBlocked)(lr->Top,nrow,irowm,ncol,icolm,y,addv);CHKERRQ(ierr); 50 IndexSpaceRestore(buf,nrow,ncol,irowm,icolm); 51 PetscFunctionReturn(0); 52 } 53 54 #undef __FUNCT__ 55 #define __FUNCT__ "MatSetValuesBlockedLocal_LocalRef_Scalar" 56 static PetscErrorCode MatSetValuesBlockedLocal_LocalRef_Scalar(Mat A,PetscInt nrow,const PetscInt irow[],PetscInt ncol,const PetscInt icol[],const PetscScalar y[],InsertMode addv) 57 { 58 Mat_LocalRef *lr = (Mat_LocalRef*)A->data; 59 PetscErrorCode ierr; 60 PetscInt bs,buf[4096],*irowm,*icolm; 61 62 PetscFunctionBegin; 63 ierr = MatGetBlockSize(A,&bs);CHKERRQ(ierr); 64 IndexSpaceGet(buf,nrow*bs,ncol*bs,irowm,icolm); 65 BlockIndicesExpand(nrow,irow,bs,irowm); 66 BlockIndicesExpand(ncol,icol,bs,icolm); 67 ierr = ISLocalToGlobalMappingApply(A->rmap->mapping,nrow*bs,irowm,irowm);CHKERRQ(ierr); 68 ierr = ISLocalToGlobalMappingApply(A->cmap->mapping,ncol*bs,icolm,icolm);CHKERRQ(ierr); 69 ierr = (*lr->SetValues)(lr->Top,nrow*bs,irowm,ncol*bs,icolm,y,addv);CHKERRQ(ierr); 70 IndexSpaceRestore(buf,nrow*bs,ncol*bs,irowm,icolm); 71 PetscFunctionReturn(0); 72 } 73 74 #undef __FUNCT__ 75 #define __FUNCT__ "MatSetValuesLocal_LocalRef_Scalar" 76 static PetscErrorCode MatSetValuesLocal_LocalRef_Scalar(Mat A,PetscInt nrow,const PetscInt irow[],PetscInt ncol,const PetscInt icol[],const PetscScalar y[],InsertMode addv) 77 { 78 Mat_LocalRef *lr = (Mat_LocalRef*)A->data; 79 PetscErrorCode ierr; 80 PetscInt buf[4096],*irowm,*icolm; 81 82 PetscFunctionBegin; 83 IndexSpaceGet(buf,nrow,ncol,irowm,icolm); 84 ierr = ISLocalToGlobalMappingApply(A->rmap->mapping,nrow,irow,irowm);CHKERRQ(ierr); 85 ierr = ISLocalToGlobalMappingApply(A->cmap->mapping,ncol,icol,icolm);CHKERRQ(ierr); 86 ierr = (*lr->SetValues)(lr->Top,nrow,irowm,ncol,icolm,y,addv);CHKERRQ(ierr); 87 IndexSpaceRestore(buf,nrow,ncol,irowm,icolm); 88 PetscFunctionReturn(0); 89 } 90 91 #undef __FUNCT__ 92 #define __FUNCT__ "ISL2GCompose" 93 /* Compose an IS with an ISLocalToGlobalMapping to map from IS source indices to global indices */ 94 static PetscErrorCode ISL2GCompose(IS is,ISLocalToGlobalMapping ltog,ISLocalToGlobalMapping *cltog) 95 { 96 PetscErrorCode ierr; 97 const PetscInt *idx; 98 PetscInt m,*idxm; 99 100 PetscFunctionBegin; 101 PetscValidHeaderSpecific(is,IS_CLASSID,1); 102 PetscValidHeaderSpecific(ltog,IS_LTOGM_CLASSID,2); 103 PetscValidPointer(cltog,3); 104 ierr = ISGetLocalSize(is,&m);CHKERRQ(ierr); 105 ierr = ISGetIndices(is,&idx);CHKERRQ(ierr); 106 ierr = PetscMalloc1(m,&idxm);CHKERRQ(ierr); 107 if (ltog) { 108 ierr = ISLocalToGlobalMappingApply(ltog,m,idx,idxm);CHKERRQ(ierr); 109 } else { 110 ierr = PetscMemcpy(idxm,idx,m*sizeof(PetscInt));CHKERRQ(ierr); 111 } 112 ierr = ISLocalToGlobalMappingCreate(PetscObjectComm((PetscObject)is),1,m,idxm,PETSC_OWN_POINTER,cltog);CHKERRQ(ierr); 113 ierr = ISRestoreIndices(is,&idx);CHKERRQ(ierr); 114 PetscFunctionReturn(0); 115 } 116 117 #undef __FUNCT__ 118 #define __FUNCT__ "ISL2GComposeBlock" 119 static PetscErrorCode ISL2GComposeBlock(IS is,ISLocalToGlobalMapping ltog,ISLocalToGlobalMapping *cltog) 120 { 121 PetscErrorCode ierr; 122 const PetscInt *idx; 123 PetscInt m,*idxm,bs; 124 125 PetscFunctionBegin; 126 PetscValidHeaderSpecific(is,IS_CLASSID,1); 127 PetscValidHeaderSpecific(ltog,IS_LTOGM_CLASSID,2); 128 PetscValidPointer(cltog,3); 129 ierr = ISBlockGetLocalSize(is,&m);CHKERRQ(ierr); 130 ierr = ISBlockGetIndices(is,&idx);CHKERRQ(ierr); 131 ierr = ISLocalToGlobalMappingGetBlockSize(ltog,&bs);CHKERRQ(ierr); 132 ierr = PetscMalloc1(m,&idxm);CHKERRQ(ierr); 133 if (ltog) { 134 ierr = ISLocalToGlobalMappingApplyBlock(ltog,m,idx,idxm);CHKERRQ(ierr); 135 } else { 136 ierr = PetscMemcpy(idxm,idx,m*sizeof(PetscInt));CHKERRQ(ierr); 137 } 138 ierr = ISLocalToGlobalMappingCreate(PetscObjectComm((PetscObject)is),bs,m,idxm,PETSC_OWN_POINTER,cltog);CHKERRQ(ierr); 139 ierr = ISBlockRestoreIndices(is,&idx);CHKERRQ(ierr); 140 PetscFunctionReturn(0); 141 } 142 143 #undef __FUNCT__ 144 #define __FUNCT__ "MatDestroy_LocalRef" 145 static PetscErrorCode MatDestroy_LocalRef(Mat B) 146 { 147 PetscErrorCode ierr; 148 149 PetscFunctionBegin; 150 ierr = PetscFree(B->data);CHKERRQ(ierr); 151 PetscFunctionReturn(0); 152 } 153 154 155 #undef __FUNCT__ 156 #define __FUNCT__ "MatCreateLocalRef" 157 /*@ 158 MatCreateLocalRef - Gets a logical reference to a local submatrix, for use in assembly 159 160 Not Collective 161 162 Input Arguments: 163 + A - Full matrix, generally parallel 164 . isrow - Local index set for the rows 165 - iscol - Local index set for the columns 166 167 Output Arguments: 168 . newmat - New serial Mat 169 170 Level: developer 171 172 Notes: 173 Most will use MatGetLocalSubMatrix() which returns a real matrix corresponding to the local 174 block if it available, such as with matrix formats that store these blocks separately. 175 176 The new matrix forwards MatSetValuesLocal() and MatSetValuesBlockedLocal() to the global system. 177 In general, it does not define MatMult() or any other functions. Local submatrices can be nested. 178 179 .seealso: MatSetValuesLocal(), MatSetValuesBlockedLocal(), MatGetLocalSubMatrix(), MatCreateSubMatrix() 180 @*/ 181 PetscErrorCode MatCreateLocalRef(Mat A,IS isrow,IS iscol,Mat *newmat) 182 { 183 PetscErrorCode ierr; 184 Mat_LocalRef *lr; 185 Mat B; 186 PetscInt m,n; 187 PetscBool islr; 188 189 PetscFunctionBegin; 190 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 191 PetscValidHeaderSpecific(isrow,IS_CLASSID,2); 192 PetscValidHeaderSpecific(iscol,IS_CLASSID,3); 193 PetscValidPointer(newmat,4); 194 *newmat = 0; 195 196 ierr = MatCreate(PETSC_COMM_SELF,&B);CHKERRQ(ierr); 197 ierr = ISGetLocalSize(isrow,&m);CHKERRQ(ierr); 198 ierr = ISGetLocalSize(iscol,&n);CHKERRQ(ierr); 199 ierr = MatSetSizes(B,m,n,m,n);CHKERRQ(ierr); 200 ierr = PetscObjectChangeTypeName((PetscObject)B,MATLOCALREF);CHKERRQ(ierr); 201 ierr = MatSetUp(B);CHKERRQ(ierr); 202 203 B->ops->destroy = MatDestroy_LocalRef; 204 205 ierr = PetscNewLog(B,&lr);CHKERRQ(ierr); 206 B->data = (void*)lr; 207 208 ierr = PetscObjectTypeCompare((PetscObject)A,MATLOCALREF,&islr);CHKERRQ(ierr); 209 if (islr) { 210 Mat_LocalRef *alr = (Mat_LocalRef*)A->data; 211 lr->Top = alr->Top; 212 } else { 213 /* This does not increase the reference count because MatLocalRef is not allowed to live longer than its parent */ 214 lr->Top = A; 215 } 216 { 217 ISLocalToGlobalMapping rltog,cltog; 218 PetscInt abs,rbs,cbs; 219 220 /* We will translate directly to global indices for the top level */ 221 lr->SetValues = MatSetValues; 222 lr->SetValuesBlocked = MatSetValuesBlocked; 223 224 B->ops->setvalueslocal = MatSetValuesLocal_LocalRef_Scalar; 225 226 ierr = ISL2GCompose(isrow,A->rmap->mapping,&rltog);CHKERRQ(ierr); 227 if (isrow == iscol && A->rmap->mapping == A->cmap->mapping) { 228 ierr = PetscObjectReference((PetscObject)rltog);CHKERRQ(ierr); 229 cltog = rltog; 230 } else { 231 ierr = ISL2GCompose(iscol,A->cmap->mapping,&cltog);CHKERRQ(ierr); 232 } 233 ierr = MatSetLocalToGlobalMapping(B,rltog,cltog);CHKERRQ(ierr); 234 ierr = ISLocalToGlobalMappingDestroy(&rltog);CHKERRQ(ierr); 235 ierr = ISLocalToGlobalMappingDestroy(&cltog);CHKERRQ(ierr); 236 237 ierr = MatGetBlockSize(A,&abs);CHKERRQ(ierr); 238 ierr = ISGetBlockSize(isrow,&rbs);CHKERRQ(ierr); 239 ierr = ISGetBlockSize(iscol,&cbs);CHKERRQ(ierr); 240 if (rbs == cbs) { /* submatrix has block structure, so user can insert values with blocked interface */ 241 ierr = PetscLayoutSetBlockSize(B->rmap,rbs);CHKERRQ(ierr); 242 ierr = PetscLayoutSetBlockSize(B->cmap,cbs);CHKERRQ(ierr); 243 if (abs != rbs || abs == 1) { 244 /* Top-level matrix has different block size, so we have to call its scalar insertion interface */ 245 B->ops->setvaluesblockedlocal = MatSetValuesBlockedLocal_LocalRef_Scalar; 246 } else { 247 /* Block sizes match so we can forward values to the top level using the block interface */ 248 B->ops->setvaluesblockedlocal = MatSetValuesBlockedLocal_LocalRef_Block; 249 250 ierr = ISL2GComposeBlock(isrow,A->rmap->mapping,&rltog);CHKERRQ(ierr); 251 if (isrow == iscol && A->rmap->mapping == A->cmap->mapping) { 252 ierr = PetscObjectReference((PetscObject)rltog);CHKERRQ(ierr); 253 cltog = rltog; 254 } else { 255 ierr = ISL2GComposeBlock(iscol,A->cmap->mapping,&cltog);CHKERRQ(ierr); 256 } 257 ierr = MatSetLocalToGlobalMapping(B,rltog,cltog);CHKERRQ(ierr); 258 ierr = ISLocalToGlobalMappingDestroy(&rltog);CHKERRQ(ierr); 259 ierr = ISLocalToGlobalMappingDestroy(&cltog);CHKERRQ(ierr); 260 } 261 } 262 } 263 *newmat = B; 264 PetscFunctionReturn(0); 265 } 266