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