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