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