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 = ISLocalToGlobalMappingApply(A->rmap->bmapping,nrow,irow,irowm);CHKERRQ(ierr); 48 ierr = ISLocalToGlobalMappingApply(A->cmap->bmapping,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),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; 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 = PetscMalloc1(m,&idxm);CHKERRQ(ierr); 132 if (ltog) { 133 ierr = ISLocalToGlobalMappingApply(ltog,m,idx,idxm);CHKERRQ(ierr); 134 } else { 135 ierr = PetscMemcpy(idxm,idx,m*sizeof(PetscInt));CHKERRQ(ierr); 136 } 137 ierr = ISLocalToGlobalMappingCreate(PetscObjectComm((PetscObject)is),m,idxm,PETSC_OWN_POINTER,cltog);CHKERRQ(ierr); 138 ierr = ISBlockRestoreIndices(is,&idx);CHKERRQ(ierr); 139 PetscFunctionReturn(0); 140 } 141 142 #undef __FUNCT__ 143 #define __FUNCT__ "MatDestroy_LocalRef" 144 static PetscErrorCode MatDestroy_LocalRef(Mat B) 145 { 146 PetscErrorCode ierr; 147 148 PetscFunctionBegin; 149 ierr = PetscFree(B->data);CHKERRQ(ierr); 150 PetscFunctionReturn(0); 151 } 152 153 154 #undef __FUNCT__ 155 #define __FUNCT__ "MatCreateLocalRef" 156 /*@ 157 MatCreateLocalRef - Gets a logical reference to a local submatrix, for use in assembly 158 159 Not Collective 160 161 Input Arguments: 162 + A - Full matrix, generally parallel 163 . isrow - Local index set for the rows 164 - iscol - Local index set for the columns 165 166 Output Arguments: 167 . newmat - New serial Mat 168 169 Level: developer 170 171 Notes: 172 Most will use MatGetLocalSubMatrix() which returns a real matrix corresponding to the local 173 block if it available, such as with matrix formats that store these blocks separately. 174 175 The new matrix forwards MatSetValuesLocal() and MatSetValuesBlockedLocal() to the global system. 176 In general, it does not define MatMult() or any other functions. Local submatrices can be nested. 177 178 .seealso: MatSetValuesLocal(), MatSetValuesBlockedLocal(), MatGetLocalSubMatrix(), MatCreateSubMatrix() 179 @*/ 180 PetscErrorCode MatCreateLocalRef(Mat A,IS isrow,IS iscol,Mat *newmat) 181 { 182 PetscErrorCode ierr; 183 Mat_LocalRef *lr; 184 Mat B; 185 PetscInt m,n; 186 PetscBool islr; 187 188 PetscFunctionBegin; 189 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 190 PetscValidHeaderSpecific(isrow,IS_CLASSID,2); 191 PetscValidHeaderSpecific(iscol,IS_CLASSID,3); 192 PetscValidPointer(newmat,4); 193 *newmat = 0; 194 195 ierr = MatCreate(PETSC_COMM_SELF,&B);CHKERRQ(ierr); 196 ierr = ISGetLocalSize(isrow,&m);CHKERRQ(ierr); 197 ierr = ISGetLocalSize(iscol,&n);CHKERRQ(ierr); 198 ierr = MatSetSizes(B,m,n,m,n);CHKERRQ(ierr); 199 ierr = PetscObjectChangeTypeName((PetscObject)B,MATLOCALREF);CHKERRQ(ierr); 200 ierr = MatSetUp(B);CHKERRQ(ierr); 201 202 B->ops->destroy = MatDestroy_LocalRef; 203 204 ierr = PetscNewLog(B,&lr);CHKERRQ(ierr); 205 B->data = (void*)lr; 206 207 ierr = PetscObjectTypeCompare((PetscObject)A,MATLOCALREF,&islr);CHKERRQ(ierr); 208 if (islr) { 209 Mat_LocalRef *alr = (Mat_LocalRef*)A->data; 210 lr->Top = alr->Top; 211 } else { 212 /* This does not increase the reference count because MatLocalRef is not allowed to live longer than its parent */ 213 lr->Top = A; 214 } 215 { 216 ISLocalToGlobalMapping rltog,cltog; 217 PetscInt abs,rbs,cbs; 218 219 /* We will translate directly to global indices for the top level */ 220 lr->SetValues = MatSetValues; 221 lr->SetValuesBlocked = MatSetValuesBlocked; 222 223 B->ops->setvalueslocal = MatSetValuesLocal_LocalRef_Scalar; 224 225 ierr = ISL2GCompose(isrow,A->rmap->mapping,&rltog);CHKERRQ(ierr); 226 if (isrow == iscol && A->rmap->mapping == A->cmap->mapping) { 227 ierr = PetscObjectReference((PetscObject)rltog);CHKERRQ(ierr); 228 cltog = rltog; 229 } else { 230 ierr = ISL2GCompose(iscol,A->cmap->mapping,&cltog);CHKERRQ(ierr); 231 } 232 ierr = MatSetLocalToGlobalMapping(B,rltog,cltog);CHKERRQ(ierr); 233 ierr = ISLocalToGlobalMappingDestroy(&rltog);CHKERRQ(ierr); 234 ierr = ISLocalToGlobalMappingDestroy(&cltog);CHKERRQ(ierr); 235 236 ierr = MatGetBlockSize(A,&abs);CHKERRQ(ierr); 237 ierr = ISGetBlockSize(isrow,&rbs);CHKERRQ(ierr); 238 ierr = ISGetBlockSize(iscol,&cbs);CHKERRQ(ierr); 239 if (rbs == cbs) { /* submatrix has block structure, so user can insert values with blocked interface */ 240 ierr = PetscLayoutSetBlockSize(B->rmap,rbs);CHKERRQ(ierr); 241 ierr = PetscLayoutSetBlockSize(B->cmap,cbs);CHKERRQ(ierr); 242 if (abs != rbs || abs == 1) { 243 /* Top-level matrix has different block size, so we have to call its scalar insertion interface */ 244 B->ops->setvaluesblockedlocal = MatSetValuesBlockedLocal_LocalRef_Scalar; 245 } else { 246 /* Block sizes match so we can forward values to the top level using the block interface */ 247 B->ops->setvaluesblockedlocal = MatSetValuesBlockedLocal_LocalRef_Block; 248 249 ierr = ISL2GComposeBlock(isrow,A->rmap->bmapping,&rltog);CHKERRQ(ierr); 250 if (isrow == iscol && A->rmap->bmapping == A->cmap->bmapping) { 251 ierr = PetscObjectReference((PetscObject)rltog);CHKERRQ(ierr); 252 cltog = rltog; 253 } else { 254 ierr = ISL2GComposeBlock(iscol,A->cmap->bmapping,&cltog);CHKERRQ(ierr); 255 } 256 ierr = MatSetLocalToGlobalMappingBlock(B,rltog,cltog);CHKERRQ(ierr); 257 ierr = ISLocalToGlobalMappingDestroy(&rltog);CHKERRQ(ierr); 258 ierr = ISLocalToGlobalMappingDestroy(&cltog);CHKERRQ(ierr); 259 } 260 } 261 } 262 *newmat = B; 263 PetscFunctionReturn(0); 264 } 265